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Cl-. Multispectral Scanner System Parameter Study and 

, 1 

Analysis Software System Description 


1. INTRODUCTION 

The utilization of sensors on earth orbiting platforms as the main 
element of an Earth Observational System has undergone substantial 
growth in recent years. ERTS-1 (Landsat-1) followed by Landsat-2 and 
-3 have proven exceptionally successful in collecting data to help 
monitor the Earth's resources. 

The principal data collection unit aboard the first three Landsats 
is the multispectral scanner known as MSS. Although this scanner has 
been providing data with a quality which exceeded most prelaunch 
expectations, it has been clear from the beginning that MSS does not 
represent the ultimate in multispectral instruments; more advanced 
instruments providing greater detail would be needed as the user community 
begins to become familiar with the use of such space data. 

The design of a multispectral scanner is a very complex matter; 
many different, interacting . factors must be properly taken into account. 
Currently operational systems such as MSS have been designed primarily 
using subjective judgements based upon experience with experimental 
data. In designing a scanner the use of empirical methods, at least in 
part, is essential. Each of the large collections of scenes which a 
given scanner will be used upon is a very complex information source; not 
enough is known to make a simple (or even a complex) model of it by 
which to make the design of a scanner a simple straightforward exercise 
of a mathematical procedure. 

And yet, more is known than when MSS was designed, and it is impor- 
tant to be able to carry out future designs on a more objective basis 
than in the past. Thus the purpose of the present work is the develop- 
ment of appropriate mathematical design machinery within a theoretical 
framework to allow: (a) formulation of an optimum multispectral scanner 

* The work in this report was done under Task 2.2C1 Multisensor Parametric 
Evaluation and Radiometric Correction Model. 
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system according to defined conditions of optimality and (b) an ability 
for convenient manipulation of candidate system parameters so as to 
permit comparison of the theoretically optimum designs -with that of 
practical approximations to it. 

In order to deal with the complexity of the design situation, the 
first step is to determine a suitable set of parameters which adequately 
characterize it but is not so large as to be unmanageable.- It has been 
observed [l] that there are five major categories of parameters which 
are significant to the representation of information in remotely sensed 
data. They are: 

1. The spatial sampling scheme 

2. The spectral sampling scheme 

3. The signal- to-noise ratio 

4. The ancillary data type and amount 

5. The informational classes desired 

Thus, it is necessary to have present in the design machinery, some means 
for evaluating the impact of change in parameter values in each of these 
five categories. 

Such a scanner design tool has been assembled in the form of a 
software package for a general purpose computer. Each of the parts of 
this package, called Unified Scanner Analysis Package (USAP) has been 
carefully devised and the theory related -to it fully documented [2, 3, 4, 5]. 
The goal of this report is to provide a documentation and description of 
the software. In constructing this documentation it was assumed that this 
package will be useful for sometime into the future, however it was also 
assumed that it will only be used by a small number of highly knowledge- 
able scientists. 

Section 2 recaps the theoretical concepts behind some of the primary 
components of USAP. These ar-e divided into (a) scanner spatial character- 
istics modeling and noise effects, (b) optimum spectral basis function 
calculations, (c) analytical classification accuracy predictions 
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(d) stratified posterior classification estimation and (e) an information 
theory approach to band selection. Although (e) is not a part of the 
USAP system, the results from this approach are helpful in understanding 
the scanner design problem. 

Section 3 shows the integration of the above modules into the 
software system. Section 4 is the user's guide to USAP describing the 
required inputs and the available output products. A listing of all 
programs is provided in the appendix. 

The work which led to USAP was immediately preceded by a simulation 
study of possible parametric values for the Thematic Mapper, a new 
scanner now being constructed for launch on Landsat-D in 1981. The 
purpose of this simulation was to compare the performance for several 
proposed sets of parameters. We will conclude this introductory section 
by briefly describing this work because it provides useful background 
and serves well to illustrate the problem. A more complete description 
of this simulation study is contained in [l, 6]. 

The general scheme used was to simulate the desired spaceborne 
scanner parameter sets by linearly combining pixels and bands from (higher 
resolution) airborne scanner data to form simulated pixels, adding 
noise as needed to simulate the desired S/N; the data so constructed 
was then classified using a Gaussian maximum likelihood classifier and 
the performance measured. The problem was viewed as a search of the 
five dimensional parameter space defined above with the study localized 
around the proposed Thematic Mapper parameters. The scope of the 
investigation was primarily limited to three parameters (a) spatial 
resolution, (b) noise level and (c) spectral bands. Probability of 
correct classification and per cent correct area proportion estimation 
for each class were the performance criteria used. The major conclusions 
from the study are as follows: 
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1. There was a very small but consistent increase in identification 
accuracy as the IFOV was enlarged. This is presumed to stem 
primarily from the small increase in signal- to-noise ratio with 
increasing IFOV, Figure 1. 

2. There was a more significant decrease in the mensuration accuracy 
as the IFOV was enlarged, Figure 2. 

3. The noise parameter study proved somewhat inconclusive due to 
the greater amount of noise present in the original data than 
desired. For example, viewing Figure 3 moving frop right to 
left, it is seen that the classification performance continues 
to improve as the amount of noise added is decreased until the 
point is reached where the noise added approximately equals that 
already initially present.* Thus, it is difficult to say for 
what signal- to-noise ratio a point of diminishing return would 
have been reached had the initial noise not been present. 

4. The result of the spectral band classification studies may also 
be clouded by the noise originally present in the data. The 
relative amount of that change in performance due to using 
different combinations of the .45— .52 pm, .74— .80 pm, .80-. 91 pm 
and .74-. 91 pm bands is slight but there appears to be a slight 
preference for the .45-. 52 pm band. The performance improvement 
of the Thematic Mapper channels over those approximating Landsa-t-1 
and -2 is clear however. 

5. Using spectrometer data it was verified that the .74-. 80 pm and 
.80-. 91 pm bands are highly correlated. 

6. Correlation studies also showed that the range from 1.0-1. 3 pm 

is likely to be an important area in discriminating between earth 
surface features. Further, it is noted that the absolute 
calibration procedure described above results in a global 
atmosphere correction of a linear type in that assuming a 
uniform atmosphere over the test site, the calibration 
procedure permits a digital count number at the airborne 
scanner output to be related directly to the present reflectance 
of a scene element. 


* The noise level in the original A/ C data was equivalent to about . 005 NEAp 
on the abscissa. See Reference [l ] - 
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Figure 1. Classification Performance vs. Spatial Resolution. 



Figure 2. RMS Error of Proportion Estimates vs. Spatial 
Resolution. 
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2. SCANNER PARAMETERS ANALYSIS TECHNIQUES 

Based upon the parametric approach introduced above, the development 
of a parametric scanner model must give explicit concern for the spatia] , 
spectral and noise characteristics of the systems. This is what has been 
done in the Unified Scanner Analysis Package (USAP) shown in Figure 4. 

USAP is composed of two distinct subsystems. The spatial aspect of it 
contains (a) a data spatial correlation analyzer, (b) a scanner IFOV 
model and (c) a random noise model. The spectral techniques are capable 
of producing an optimum spectral representation by modeling the scene as 
a random process as a function of wavelength, followed by the determination 
of optimum generalized spectral basis functions. Conventional spectral 
bands can also be generated. Also studied was an information theory 
approach using maximization of the mutual information between the reflected 
and received (noisy) energy. The effect of noise in the data can be 
simulated in the spectral and spatial characteristics. Two different data 
bases are used in the system. The spectral techniques require field 
spectral data while the spatial techniques require MSS generated data, 
aircraft and/or satellite. The system performance, defined in terms of 
the classification accuracy, is evaluated by two parametric algorithms. 

A detailed system description and user’s guide is presented in Sections 3 
and 4. In the following, the theoretical ideas behind the five major 
elements of USAP are discussed. 

2.1 Analytical Classification Accuracy Prediction 

Throughout the analysis of remotely sensed data, the probability of 
correct classification has ranked high among the set of performance indices 
available to the analyst. This is particularly true in a scanner system 
modeling where generally the optimization of various system parameters 
has as its prime objective the maximization of the classification accuracy 
of various classes present in the data set. 

The estimation of the classification accuracy is fairly straightfor- 
ward if Monte-Carlo type methods are employed. In system simulation and 
modeling however, such approaches are generally a handicap due to their 




Figure 4 


Block Diagram of the Unified Scanner' Analysis Package (USAP) . 
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heavy dependence on an experimental data base, the availability of which 
can be limited due to a variety of reasons. What is required, therefore, 
is a parametric classification accuracy estimator for a multiclass, 
multidimensional Gaussian Bayes classifier. This procedure should require 
the class statistics, mean vectors and covariance matrices, as its only 
input and produce a set of probabilities of correct classification. 

This technique has been developed, tested, implemented and comprehensively 
reported [2], The following is a summary of the method and some' results. 

The probability of Error as an N-Tuple Integral . 

The classification of a multidimensional observation vector into 
one of M populations is conceptually identical to the binary case. Let 
ft, M and N be the feature space, number of classes and the dimensionality 
of ft respectively. The procedure is to divide ft into M mutually disjoint 
sets, I*., and to assign each feature vector to a set in accordance with 
an appropriate rule. This is illustrated in Figure 5. Let Z_^, i •= 1, 

2, ..., M partition ft in R K . The Bayes risk is defined as 


R = 


M 

’ l 

i=l 


'Z. 


M 

l 

J=1 


POO 


-ij 


f'(xk) 


dx 


( 1 ) 


where C is the cost of deciding w. where m. is true. In the case where 
ij 13 

C =0 for i=i and C..=l for j . R is the probability of error, 
ij 13 


Among all possible choices of Z. the Bayes rule partitions ft into 

Z =Z * such that R=R* is the the minimum probability of error. Assuming 
i i 

that the population statistics follows a multivariate normal law, the 
optimum Bayes rule is as follows [7]. 


X e w. if W. < W. Vi i j = 1, 2, ..., M 
- i 13 


-1 

W. = ( X - P.) T j. (X - jj.) + ln|L |-21n P(co ) 
1 1 1 1 


where 
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with 


X = observation vector 

11 . = mean vector for class to , 

— i i 

T , = covariance matrix for class to . 
“i 1 

P(to ) = apriori probability for ok ' 


( 2 ) 


The error estimate based on direct evaluation of Eq. (1) exhibits all the 
desired, properties outlined previously. 


The evaluation of multiple integrals bears little resemblance to 
their one dimensional counterparts, mainly due to the vastly different 
domains of integration. Whereas there are three distinct regions in one 
dimension; finite, singly infinite, and double infinite; in an N dimensional 
space there can potentially be an infinite variation of domains. The 
established one dimensional integration techniques, therefore, do not 
carry over to N dimensions in general. Hence, it is not surprising that 
no systematic technique exists for the evaluation of multivariate integrals 
except for the case of special integrands and domains The major 

complicating factor is the decision boundaries defined by Eq. (2). T is 
defined by a set of intersecting hyperquadratics. Any attempt to solve 
for the coordinates of intersection and their use as the integration limits 
will be frustrated if not due to the cumbersome mathematics, because of 
impractically complicated results. 


In order to alleviate the need for the precise knowledge of boundary 
locations and reduce the dimensionality of the integral, a coordinate 
transformation followed by a feature space sampling technique is adopted. 
The purpose of the initial orthogonal transformation of the coordinates 
is an N to 1 dimensionality reduction such that the N~tuple integral is 
reduced to a product of N one dimensional integrals. Let the conditional 

classification accuracy estimate, P i be the desired quantity. Then 

' i 

the transformed class to. statistics is given by 
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M. = 0? IK 3 = 1» 2, ..., M (3) 

—3 — * 3 

S . - 0 1 Zj 0 

where 0. is the eigenvector matrix derived from Naturally, in' each 

transformed space, T^ft), ok has a null mean vector and a diagonal 
covariance matrix. 

The discrete feature space approach is capable of eliminating the 

need for the simultaneous solution of M quadratic forms. If ft is the 

continuous probability space, a transformation T_^ is required such that 

in T (ft), T can be completely described in a nonparametric form, thereby 
i i 

bypassing the requirement for -an algebraic representation of I\. This 

desired transformation would sample ft into a grid of N— dimensional 

hypercubes. Since the multispectral data is generally modeled by a 

multivariate normal random process, a discrete equivalent of normal 

random variables that would exhibit desirable limiting properties- is 

required. Let y ~Bi (n, p) be a binomial random variable with parameters 

n and p. The x defined by 
1 n 


x 


y n - np 


ll /np (1-p) 


“ 0, X , 2 , n 


(4) 


converges to x~N(0, 1) in distribution [9], i.e., 


lim F (X) *■ F(x) 

n 


The convergence is most rapid for p=h, then 


x 


n 


(y n - n/2)2 
/n 


(5) 
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The variance of x is set equal to the eigenvalues of the transformed Y. 

n 1 

by incorporating a multiplicative factor in Eq. (5). 

The segmentation of by a union of elementary hypercubes makes 
nonparametric representation of I\ and its contours feasible. Following 
the orthonormal transformation on w_. and sampling of 0 . accordingly, each 
cell's coordinate is assigned to an appropriate partition of F . This 
process is carried out exhaustively, therefore can be defined as a set 
such that 


r. 

i 


= {UX :X er.} 

— n — n i 


(6) 


once the exhaustive process of 
f(x|a. ) over r i is represented 
elementary cells within F^. The 


C. 

l 


f 


(x|oO 



(x i “i ) 


assignment is completed, the integral of 
by the sum of hypervolumes over the 
elementary unit of probability is given by 



(x 2 K ) 



P" t 

fa 

2 (7) 


where C. is the domain of a sampling cell centered at the origin and 6^ 
is the width of a cell along the ith feature axis. The conditional 
probability of correct classification is therefore given by 


P to, 
c 1 


1 cL J 


C 1 + 6 1 


2 f (x-Jok) I i (C) dx x 


c 2 + 6 2 


2 f (x 2 |to i ) I (C) dx 2 


c r 6 ! 


c 2~ 6 2 


( 8 ) 


c +6 
n n 


2 £ (k n |».) I. (C) 


c -6 
n n 


with overall classification accuracy given by 
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P 

c 


M 


i-1 


P(w.) P 
1 . c 


(9) 


where 


I i (C) 


1 if CeT. 

i 

0 otherwise. 


C = The domain of an elementary cell 


( 10 ) 


Pigure 6 is a geometrical representation of Eq. (8). 

Experimental Results . 

The analytic classification accuracy prediction (ACAP) has been 

thoroughly tested and documented [2]. Two examples are repeated here. 

The first experiment investigates the performance of the estimator vs. 

grid size i.e., number of cells per axis, n. Small to moderate range of 

n is required if computation time is to remain realistic. Figure 7 

shows the variation of P i vs. n for three classes having some hypothe- 

c w . 

' x 

tical statistics in 3 dimensions. The main property of the estimator is 
its rapid convergence toward a steady state value thereby alleviating 
the need for excessively fine grids and hence high computation costs . 


The data collected over Graham County, Kansas is used to perform 
a comparison between the ACAP algorithm and a ratio estimator such as 
LARSYS. The results are tabulated in Table 1. 




NO OF CELLS PER AXIS 


Figure 7. ACAP Glassification Accuracy Estimate vs. Grid Size for the Test Population Data 
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Table 1. ACAP and LARSYS Performance Comparison, 


Class 

No. of Points 

LARSYS 

ACAP 

Bare Soil 

443 

65.9 

78.3 

Corn/Sorghum 

99 

89.9 

91.0 

Pasture 

1376 

98.4 

95.0 

Wheat 

459 

94.8 

93.9 

Overall 

2377 

87.2 

89.6 


The comparison of ACAP and LARSYS results from Table 1 is inconclusive. 

In some cases the difference is negligible (corn) and in some, significant 
(bare soil) . Examination of the data statistics revealed that the 
assumption of normality does not hold throughout the populations statistics. 
This problem can be rectified by simulating random Gaussian data having 
identical statistics with the real data, hence assuring the normality 
assumption. Repeating the LARSYS and ACAP procedures produces a new set 
of classification accuracies, Table 2. 

Table 2. ACAP and LARSYS Performance Comparison Simulated Data. 


Class 

LARSYS % 

ACAP% 

Accuracy Difference% 

Bare Soil 

77.8 

78.3 

0.5 

Com 

91.2 

91.0 • 

0.2 

Pasture 

95.3 

95.1 

0.2 

Wheat 

94.2 

93.9 

0.2 

Overall 

89.6 

89.6 

0 


The differential between ACAP and LARSYS results has been narrowed 
considerably, ranging from a high of 0.5% for bare soil to 0% for the 
overall classification accuracy. Tito conclusions can be drawn from the 
results of this experiment. First, the ACAP and Monte Carlo type 
classifiers produce practically identical results if the underlying 
assumptions are satisfied (e.g., normality of the statistics). Second 
and more revealing is the fact that the results of the ACAP processor 
indicated an upper bound for the classifability of bare soil had 
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its statistics been Gaussian. This result is a direct property of ACAP's 
data independence. Figure 8 is the ACAP estimator vs. n for Graham Co. 
data. 


The above selected experiments and others reported in the biblio- 
graphy establishes ACAP as a viable and necessary tool in any analytical 
remote sensing data collection system modeling and simulation when the 
performance index is defined as the probability of correct classification. 

2.2 Stratified Posterior Classification Performance Estimator. 


The second classification accuracy estimator to be presented here 
(SPEST) is based on the maximization of the aposteriori probability 

l 

associated with each sample. This formulation is closely related to 
the maximum likelihood principle used in the ACAP . The distinction 
arises in the determination of integration domains. Where in ACAP a 
"deterministic" grid was set up to sample the feature space, SPEST uses 
an internally generated random data base and assigns the feature vector 
to the appropriate class via the maximum aposteriori principle. Due to 
the different approaches adopted, the statistical properties of the 
estimators could be substantially different although no major study has 
yet been carried out. It has been observed however, that the SPEST 
algorithm is somextfhat faster than the ACAP in selected cases. The 
approach here is similar to that described in Moore, Whitsitt, and 
Landgrebe [10] . 

Let X be an observation from one of M classes on, i = 1, 2, 3, . . . , M, 
with a priori probabilities P^. The maximum likelihood decision rule can 
be stated as follows: Assign X to the class if 

P(cjj.. |X) = max {P(w.|x)} 

K .1 

1 

This rule partitions the observation space 52 into subregions 
corresponding to the classes •• •» respectively. Define the 

indicator function as 
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Figure 8 


ACAP Classification Accuracy Estimate vs 


Grid Size for Graham Co., Kansas Bata. 
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I ± (x) 


1 xsQ. 
1 

0 

x 


The probability of correct classification integral is given by 


P 

c 


f M 

l P. I.(x) p(x[w ) dx 
i=l 1 1 1 


(ID 


It is desirable to evaluate the probability of correct classification for 

3 

each class as well as the overall probability. The performance probability 
for the ith class is 


P ci = J fl V X) P (aE l“j? 


The overall performance, then, is 


(12) 


P 

c 


M 

I 

i=l 


P.P . 
1 ci 


( 13 ) 


From Bayes' rule 


P(to i |x) p(x) 

P(x|u.) = ' ' p 

i 


hence. 


P . 
ci 


P(«.|x) 

I.-(x) — p(x) dx 

J Q 3 ? i 


M 

-p(x) = l P p(x|m ) 

j=l J J 


p(x) is a mixture density 
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Therefore, 


M P. r 

P . = £ I.OO P(«> |x) p(x|(D ) dx (14) 

C ‘ l j=l i .(2 1 


Define 


Q ± (x) = I 1 (x) P(to ± |x) 


Then 


j Q.(x) p(x|o) ) dx 
1 2 


is the conditional expected value of Q^(x) given that X comes from the 

class C.. The estimate 
J 


P . 
ci 


P. , “3 

. , »: J, W 

3=1 l 3 k=l 


M 

l 


(15) 


is unbiased. This estimator is similar to the stratified posterior 
estimator described by Whitsitt [lO], 

To do this a pseudo-random sequence of uniformly distributed random 
digits is generated by the power-residue method and is transformed by the 
inverse cumulative-distribution— function method to obtain nearly Gaussian 
samples. These samples are used to fill the elements of the data vector Y_. 
Each vector Y, then, has expected value _0 and covariance matrix I_. 

By performing the transformation 

k 

X = T. Y + m. 

- 3 “J — ~1 


( 16 ) 
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on the vectors Y, the random vectors for class j are obtained, where 

“3 

is the matrix of eigenvectors required to diagonalize the covariance matrix 
of class j, T is the diagonal matrix of eigenvalues and nu is the desired 
mean vector. These random vectors are used to evaluate the estimators in 
Eqs. (15) and (13). 


The term that must be evaluated from Eq. (14) is 


P p(x|o) ) 

P(C. |x) = - ■ 

I \ p(x| "w 


To evaluate this probability compute 
the largest value of the product P. 


P. p(x|oj.) for each class m.. 
3 1 J 3 

p(x w.) and divide by the sum 
J 


I 

k 


Choose 

p k p(x K } 


p(x|w ) = 


(2,r) N/2 [K k r 


— exp {-%(x-B k ) T K,. 1 (x - m^) ) 


(17) 


m^ and are the mean vector and covariance matrix respectively for class 
k. Substituting Eq* (16) into (17), 


p (x co^) = 


T „-l 


(^) W/2 |K k l 


- exp {-h [Y X r l! K X $ r 2 Y + 
v, J 3 K J J 


2X t I? d 4 1 ^ 


T „-l 


(18) 


In this form it is not necessary to perform the intermediate, computational 
step of transforming the data. We need only to generate M sets of random 
vectors Y with mean vector 0 and covariance matrix JE and use them in the 
Eq. (18). 
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Estimator Evaluation . 

A subroutine program was written to evaluate classification perfor- 
mance by the above method. To test the method a three class problem was 
constructed. The mean vectors for the classes were 

Mj_ = [-1* -1, -1] T 

M 2 - [o, 0, of 

M 3 = [1, 1, .... if 

The covariance for each class was the identity matrix. The number of 
random vectors generated for each class was 1000. The exact classifica- 
tion accuracy as a function of the dimensionality can be evaluated for 
this case 

P _ = 1 - erfc (^N/2) 
cl 

P = 1 - 2 erfc (v¥/ 2) 
cz 

P _ = 1 - erfc (v^J/2) 
c3 

P =1-4/3 erfc (^/2) 
c 

f OT £ -x 2 /2 

where erfc (a) = dx 

^ a / 2 tt 

and n is the number of dimensions. Table 3 contains the results of 
evaluating the class conditional performance and overall performance 
from one to ten dimensions. 

To evaluate the variance of the estimates different starting values 
for the random number generator were used. Twenty trials were used to 
evaluate the maximum bound and the standard deviation from the true value. 
These results are presented in Table 4. 

For the overall accuracy the estimate is within .005 of the true 
value. This is certainly sufficient for performance estimation. The 



Table 3. Test of the SPEST Error Estimate. 



P 

C 1 

P 

C 2 

P 

c 3 

P. 

c 

A 

P cl 

A 

P c2 

hcJ- > 
O 
U> 

A 

P, 

C 

1 

0.6915 

0.3829 

0.69'15 

0.5886 

0.6859 

0.3793 

0.7001 

0.5884 

2 

0.7602 

0.5205 

0.7602 

0.6803 

0.7671 

0.5116 

0.7700 

0.6829 

3 

0.8068 

0.6135 

0.8068 

0.7423 

0.8037 

0.6202 

0.8081 

0.7440 

4. 

0.8413 

0.6827 

0.8413 

0.7885 

0.8283 

0.6852 

0.8550 

0.7895 

5' 

0.8682 

O'. 7364 

O'. 8682 

0.8243 

0.8642 

6.7425 

0.8703 

0.8256 

6 

0.8897 

0.7793 

0.8897 

0.8529 

0.8767 

6.7939 

0.8787 

0.8498 

7 

0.9071 

6.8141 

6.9071 

0.8761 

6.8993 

6.8242 

0.9065 

0.8766 

8 

0.921-4 

0.8427 

0.9214 

0.8951 

0.9129 

0.8472 

0.9240 

0.8947 

9 

0.9332 

0.-8664 

6.9332 

0.9109 

0.9193 

6.8809 

0.9360 

0.-9120 

10 

0.9431 

0.8862 

6.9431 

0.9241 

0.9209 

0.90i2 

0.9481 

0.9234 


N> 

■P* 
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Table 4. Variation of Estimates. 


P 


cl 




P 

c 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


.016 

.033 

.010 
' .019 

.017 

.049 

.003 . 
.005 

.018 

.036 

.010 

.018 

.014 

.027 

.002 

.005 

.016 

.046 

.017 

.031 

.017 

.055 

.003 

.007 

.011 

.025 

.016 

.029 

.015 

.029 

.003 

.005 

.015 

.031 

.014 

.033 

.012 

.026 

.002 

.004 

.014 

.026 

.014 

.023 

.010 

.022 

.003 

.006 

.009 

.027 

.016 

.033 

-.012 

.027 

.003 

.005 

.013 

.025 

.013 

.036 

.012 

.023 

.003 

.006 

.013 

.026 

.014 

.031 

.012 

.021 

.002 

.004 

.009 

.016 

.012 

.024 

.009 

.019 

.002 

.005 


a 

Bound 

a 

Bound 

a 

Bound 

a 

Bound 

0 

Bound 

a 

Bound 

a 

Bound 

0 

Bound 

a 

Bound 

0 

Bound 


0 = standard deviation 

Bound - maximum difference between estimate and true value over 
20 trials 
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class conditional estimates are less reliable but are sufficient to 

observe trend in the performance due to individual classes. 

« * 
i 

2.3 Scanner Spatial Characteristics Modeling 

The multispectral scanner represents the. most important element in 
a remote sensing data gathering system. Therefore, an understanding of 
the signal flow through this subsystem is essential. As data is processed 
through the scanner, its statistical properties undergo a transformation. 
This in turn will alter the population separabilities and hence the 
classification accuracies. The comparison of this quantity at the scanner 
input and output and observation -of its variation with the system para- 
meters. sheds considerable light on the overall system design. Since the 
Bayes spectral classifier depends solely on the population of spectral 
statistics, methods need to be developed that relate the scanner’s input and 
output statistics. A complete derivation of such relationship is given 
in Appendix A of [ 2 3. A summary follows: 


Scanner Characteristic Function. 


s 

Figure 9 is a basic block diagram of the scanner spatial model where 
through are N stochastic processes corresponding to N .spectral 
bands and h(x,y) is a two dimensional PSF. In particular where the Landsat 
scanner' is concerned, the assumption of a Gaussian shaped IFOV has been 
widespread. Let f(x,y), g(x,y) and h(x,y) denote the input and output 
random processes associated with any two matching bands and' the scanner 
PSF respectively. It is well known that the above quantities are related 
by a convolution integral. 


g(x,y) 


* 

f(x-X^,y-A 9 ) d^dA^ 


(19)- 


it follows that 


S wo (u,v) - S, f f {u,v) jH(u,v) | 2 


( 20 ) 
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where S(u,v) is the spectral density of ;the appropriate random process 
and H(u,v) is the two dimensional Fourier transform of the scanner PSF-. 

Let M(u,.v) = |h(u, v)| 2 , and m(x,n) its inverse transform. Then the output, 
spatial autocorrelation function is- given by 


R (r,n) = ^.(Tjn) * m(t, n) (-21) 

gg ff 


In order to obtain specific results, the following assumptions are 
invoked; (a) exponential data spatial correlation, (b) Gaussian IFOV, 


R f£ (T,n) 




x,n o, 1, 2', ... 


2 .2 

, . . -x -y 

h(x,y) = c^ e •*— j e — j 

r r 

o o 


( 22 ) 


*“3, “b- 

where p = e and p = e are the adjacent sample and line correlation 
x • K y 

coefficients respectively, r^ is the scanner PSF characteristic length 
in- pixels and a constant providing unity filter gain. Using the 
separability property of the functions involved-i 


R 

gg 


(r„n) 


R^Cx-x) m(x). dx 


R f f (try) m(y) dy 


where 


2 ' 2 2 2 
c r _ x __ n 

m(x,n) = ~ — 2 ^— e 2r 2 6 2r 2 

o o 


( 23 ) 


carrying out the integration, the scanner characteristics function is 
given by 



29 


W (x,n,a,b) = 
s 


2 2 
a r 

o 


-ax 


2 2 
a r 

o 

Q(ar + e 2 

o r 

o 


+ax 


Q(ar + — ) 
or 

o 


X 


,2 2 

b r 

o 


-bn 


. 2 2 
b r 

o 

Q(br ) + e 2 
or 

o 


+bn 


Q(br + 12 -) 
^ or 

o 


(24) 


where 


Q(x) 



^ v^2tt 


The spectral statistics is a subset of the spatial statistics hence 

W (0,0,a,b) defines the ratio of the variance at the scanner output to 
s 

the corresponding input quantity. 


The output crosscovariance terms can he similarly derived. Let 
the crosscorrelation function between bands i and j be defined as 


-a M ~b In 

R f (x,n) - r f £ a £ a f> e 
x J x 3 i .1 


(25) 


where r ^ is the spectral crosscorrelation- coefficient at the input 
£ i j 

such that I r , i < 1 * a . . and b . . are defined similar to a and b with 
1 f j 1 - xj iJ 

the additional channel specification. Following the previous technique 
it follows that the crosscovariance term between channels i and j is 
given by 


Ys. <0 - 0) ■ w s (0 ’°> i! i3’ b ij ) 


. .x. . x. . x. . 1 , 1 . 
X j x j i j 


(26) 
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The corresponding crosscorelation coefficient follows : 


' § i g 3 


( 27 ) 


Therefore, the band-to-band correlation coefficients are identical at 
the scanner input and output provided the spatial auto and crosscorrela- 
tion function at the input are equivalent, i.e., a ±i =a , b^b^ • A 
closed expression can also be obtained for a rectangular PSF defined by 


h(x,y) 



M < r o/2 


[0 elsewhere 


(28.) 


the corresponding characteristic function is given by 


W (0,°,a Kp -77^-d- 

11 o 


1-e 


-a. . r 
11 o 


a. .r 
xi o 


-) X 


2 


b . .r 
11 o 


(1 


1-e 


-b . .r 
11 o 


b. .r 
11 o 


(29) 


Eq. (24) and Eq. (29) are plotted in Figures 10 through 13 for different 

scene correlations with the IFOV = t q as a running parameter. The 

universal property exhibited by W of either PSF is the increasing output 

s 

variance reduction as IFOV is increased. This property has, been widely 
verified experimentally. Comparison of Figures 11 and 13 indicates that 
for the same IFOV, a Gaussian shaped IFOV causes a heavier variance 
reduction in a spectral band than that of a rectangular PSF. This 
property can potentially produce a higher separability among the popula- 
tions as the signal is processed through the scanner electronics. 
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Figure 13. Scanner Characteristic Function vs. Scene Correlation. 
Adjacent Line Correlation = .8. 
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The input-output statistical relationship just developed along with 
the analytic classification accuracy predictor, provide the two basic 
tools required for a parametric evaluation of the MSS performance under 
varying operating conditions. As an illustration, three hypothetical 
classes with some prescribed statistics are specified at the scanner 
input. Three different sceneswith adjacent sample correlations of 0.5, 
0.85 and 0.95 are considered. The scanner characteristic function 
produced a set of transformed statistics at the output followed by the 
estimation of the conditional classification accuracies using ACAP , for 
8 different IFOV's. The results are plotted in Figures 14 through 16. 

Two main properties stand out. First is the improvement in class 
separability as the XFOV is enlarged. This is true in all the cases. 

The rate of improvement however, is strongly correlated with the scene 
spatial correlation. The lower the input scene' correlation, the higher 
the classification accuracy improvement per IFOV step. This property 
is brought about by the characteristics of W g where one step increase 
in IFOV size produces a greater variance reduction for a low scene 
spatial correlation than a similar increase would cause in a highly 
correlated scene. 

The scene spatial correlation plays a significant role in the overall 
system performance which is not readily obvious. One of the well 
known properties of linear systems with random inputs is the reduction 
of the output variance/ input variance ratio as the point spread function 
(PSF) is widened. In this section it has been shown theoretically that a 
third factor in this reduction is the spatial correlation structure of 
the input stochastic process. Specifically, with everything else fixed, 
a process having a moderate scene correlation will undergo a tighter 
clustering around its mean than an otherwise identical process with a 
highly correlated spatial characteristic. On the extreme side of the 
correlation scale with a small pixel- to-pixel correlation, the ratio of 
the output variance to the corresponding input quantity is very negligible. 
Consider a bandlimited white noise process with a spectral density, 
shown in Figure 17, where W is the bandwidth and No/2 the two sided spec- 
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SCANNER IFOV IN HIGH RESOLUTION PIXELS 


Figure 14. Scanner Output Classification Accuracy vs. I^V. 
Adiacent Sample Correlation = .55. 


LO 



SCANNER IFOV IN HIGH RESOLUTION PIXELS 


figure 15. Scanner Output Classification Accuracy vs. IFO'L 
Adiacent Sample Correlation = .8. 





SCANNER IFOV IN H I GH~ RESOLUTION PIXELS 


Figure 16. Scanner Output Classification Accuracy vs. IFOV. 
Adjacent Sample Correlation =\95. 
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No/ 2 

— 2irW 2ttW (0 

Figure 17. Bandlimited White Noise Spectral Density. 

tral density. As W increases the adjacent pixel correlation in the 
scene decreases. The increase in W, however, is accompanied by a 
decrease in No/2 if this process is to remain physically realizable 
(finite energy). Under a finite energy constraint, therefore, as W s- «, 

No/2 0. In the limit the energy content of the output random process 

will be nil. 

Random Noise . 

Additive random noise entered at various stages of a scanner system 
can degrade the overall system performance substantially. The noise can 
be classified into two broad categories! external and internal. A 
major source of external noise is atmospheric in nature mainly due to 
absorption (e.g., x-rater vapor) and scattering. The detector and quanti- 
zation noise comprise the major component of the internal noise sources. 

From a system analysis point of view, the latter represents a more tract- 
able and better understood component of the random noise [ll], while 
the former still awaits further experimental documentation. The purpose 
of this work is not so much the exploration of the various noise sources 
but the integration of its effect within an analytic analysis package 
once its characteristics and origin has been determined. 

From the theoretical results obtained it can be stated that atmospheric 
noise, in the uplink path at least, has negligible degrading factor 
compared with the detector and quantization noise. Let f(x,y)-, N^(x,y), 

£' (x,y) and N^Cxjy) be the input random process, input additive white 
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noise, the output random process and the noise component of the output 
signal respectively, then 

f"(x,y) = f(x,y) * h(x,y) (30) 

N .(x,y) = N (x,y) * h(x,y) (31) 

f 

Define 

( SNR) £ = Var {f(x,y)}/Var {N f (x,y)> (32) 

(SNR) „= Var {f " (x,y) >/Var (N „(x,y)} (33) 

f f 

Recalling the functional dependence of W g on the input scene spatial 
correlation, it follows that the ratio of the variance of a white noise 
process at the scanner output to the corresponding input quantity is of 
the order of 5% to 10%, higher or lower depending on the IFOV size. 
Therefore, 


Var {f"(x,y)l < Var (f(x,.y)'} (34) 

Var (N „(x,y)} « Var {'N f (x,y’>} (35) 


hence 


(SNR) _ » (SNR) (36) 

f 

It then follows that the noise component of the output process prior to 
detector and quantization noise is negligible in most cases. 

In order to observe the effect of noise on the scanner output class 
separability the test class statistics were modified to exhibit the effect 
of random noise. The assumed properties of the noise are additive, 
whine and Gaussian. Let r'(x,v) be the signal to be telemetered to 
Eart h . 
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f""(x,y) = f'(x,y) + N _(x,y) (37) 

f 

the statistics of f^^y) and f'(x,y) are related by 

l = l+l (38) 

-f" ~f' ' 

f 

the simple addition is due to the signal and noise independence. Assuming 

a zero mean N the mean vector are identical, i.e., 
f 

E if"} = E {f'J 

Among the four assumptions about the noise, tits Gaussian property is the 
weak link due to the Poisson distributed detector noise and uniformly 
distributed quantization noise. Relaxing the Gaussian noise assumption, 
however, would mean the design of an optimum classifier for non-normal 
classes and evaluation of its performance. A task that would complicate 
matters considerably. Due to the relatively insufficient documentation 
of the characteristics of random noise in multispectral’ data, the initial 
Gaussian assumption is adhered to. 

Following the adopted SNR definition, three different noise levels 
are considered and the corresponding overall classification accuracies 
for the three previously used test classes are estimated. Figure 18 is 
the variation of P vs. IFOV with SNR as the running parameter. For a 

A A 

fixed IFOV, P^ increases with increasing SNR. For a fixed SNR, P^ 
increases with increasing IFOV size. These illustrations have shown 
that with a proper coupling between the ACAP and the scanner characteristic 
function, the progress of the population statistics through the system 
can be studied on an analytical and entirely parametric basis. The 
accompanying classification accuracies can measure the designer's success 
in selecting the spatial and/or spectral characteristics of a Multispectral 
Scanner System. 



SCANNER IFOV IN HIGH RESOLUTION PIXELS 


Figure 18. Overall Output Classification Accuracy Variation with Noise and IFOV. 
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2.4 Optimum Spectral Function Research 

In earth observational remote sensing much work has been done with 
extracting information from the spectral variations in the electromagnetic 
energy incident on the sensor. Of primary importance for a multispectral 
sensor design is the specification of the spectral channels which sample 
the electromagnetic spectrum. An analytical technique is developed for 
designing a sensor which will be optimum for any well-defined remote 
sensing problem and against which candidate sensor systems may be compared. 

Let the surface of the earth at a given time be divided in strata 
where each stratum is defined to be the largest region which can be 
classified by a single training of the classifier. Each point in the 
stratum is mapped into a spectral response function X(A) as in Figure 19. 
That is if one observes a point in the stratum with the sensor, then the 
function X(A) describes the response variations with respect to the 
wavelength, A. The stratum together with its probabilistic description 
defines a random process, and the collection of all of such functions 
X(A) which may occur in the stratum is called an ensemble. 

The general concept of a pattern recognition system in this applica- 
tion requires that if each X(A) is to be classified by a classification 
algorithm, this can be accomplished by first measuring a finite number 
of attributes of X(A) , called features. This is the function of the 
sensor system as depicted in the upper left portion of Figure 20 where 
X^, X X^ are the values of N features for a given X(A) . It 
may be viewed as a filtering operation on X(A)-. 

For example, on the right portion of Figure 20 the function of MSS 
of the Landsat satellites is illustrated. In this case a number propor- 
tional to the average energy in a wavelength interval is reported out by 
the sensor for each of four wavelengths. Mathematically this may be 
expressed as 


X 


n 



(A)dA 


n = 1, 2, 3, 4 
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Next we must consider what would constitute an optimum sensor. We 
first note that in general the sensor may be used over any part of the 
earthfe surface, at anytime, and for many different applications (sets 
of classes). Therefore the sensor must be optimized with respect to the 
entire set of strata represented by these cases. As a result of the large 
size of this set and the fact that its statistical description is not 
known, we will optimize the sensor with respect to its signal representa- 
tion characteristics. The (X(A) } each contain information useful to the 
classifier; we require of the sensor design that for a given N a maximum 

of this information which was in X(A) still be present in (X }. Since 

n 

the specific nature of this information is not known a priori, we can 

only assure that this will be the case for any stratum if X(A) is 

recoverable from (X } . 

n 

Let X(A) be the result of attempting to reconstruct X(A) from (X }. 

n 

A fidelity criterion which is useful in this instance is 

e = j [X(A> - X(A)] 2 dA (39) 

i 

the so called mean square error or mean square difference between X(A) 
and X(A) . 

It is known [12] that a reconstruction scheme which minimizes e for 
a given N is 


X(A) = X 1 <j) 1 (A) + X 2 4 2 (\) +.... -P X N <j> n (A) 


N 


= I 

n=l 


X <t> (A) 
n n 


(40) 


provided that the {^(A)} are orthogonal over the wavelength interval of 
interest, i.e., 


<f> (A) cf> (A) dA = 0 
m n 


m f n 


(41) 




x> 
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and the {X } are calculated by 
n 


X 

n 


X(A) 4> (A) dA 

J n 


( 42 ) 


Note for example that the Landsat example of Figure 20 satisfies 
these conditions. In the lower right of Figure 20 is depicted the 
result of such a reconstruction for the Landsat example. 


While use of Eq. (42) in the case does minimize e with respect to 

the choice of values of {X }, a further improvement may by obtained by 

n 

choosing a set of { (f> . { A | } which minimizes e. It can be shown [5, 12] 

that the set { (p (A)} which accomplishes this must satisfy the equation 
n 


a .<(> . (A) 
xx 


R(A,5) * ± (C) d£ 


(43) 


where 


R(A,C) = E([X(A) - m(A)] [X(5) - m(5)]l (44) 

is the correlation function of the random process and m(A) is its mean 
value at A . 

Such a signal representation defined by Eqs . (40-44) is known as a 
Karhunen-Lohve expansion [l3]. It provides not only for the most rapid 

A 

convergence of X(A) to X(A) with respect to N but in addition the random 

variables {X } are uncorrelated and since the random process is Gaussian 
n 

they are statistically independent. Further the only statistic required 
of the ensemble is R(A,£). This representation of {X(A) } is therefore 
not only optimal, it is convenient. 

A useful generalization of the Karhunen-Loeve expansion can be made. 
Suppose a priori information concerning portions of the spectral interval 
are known and it is desired to incorporate- this knowledge into the analysis . 
A weighting function w(A) , is introduced which weights portions of the 
interval according to the a priori information. As an example, measure- 
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ments were taken over the spectrum and it was observed that there was 
considerable variation in the signal in the water absorption bands 
around 1.4 and 1.8 micrometers. This variation was due to measurement 
and calibration difficulties rather than being a result of variations 
in the scene. Therefore, the weighting function was set to zero in 
these absorption bands. This generalization is referred to as the 
weighted Karhunen-Loeve expansion [5]. Eqs. (40), (42) and (43) become 


X(A) = l x. <}> (A) 

. i w . 

1=1 i 


AeA 


(45) 


a 4> (A) = 

w. w. 
i x 


R(A,£) w(£) <J> (£) d£ 

w . 

l 


(46) 


x. = I X(A) w.(A) <J> (A) dA 

1 i w . 

X 


(47) 


where the eigenfunction, <j> (A) are solutions to the integral Eq. (46) 

* Wi. 

X 

with the weight w(A). The special case where w(A) = 1.0 for all AeA 
reduces the expansion to the original form' in Eq. (40), (42) (43), and 
(44). i 

The results of having utilized this means of optimal basis function 
scheme on spectral data are contained in reference [5]. From them one car 
see the significant improvement in classification accuracy which decreased 
spectral representation error will provide. One can also determine the 
spectral resolution and band placement needed to achieve such classification 
accuracy improvement. 

2.5 Information Theory Approach to Band Selection 

The problem of selecting a set of "optimum"' windows in the electro- 
magnetic spectrum for observing the reflected sunlight has always been of 
considerable interest. Depending on the definition of the optimality 
different methods have been developed. One such approach was shown in 
Section 2.4 using K—L expansion to select an optimum set of basis functions. 
In this section an information theoretic definition of optimality developed 
in [3] is explored. 
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Mutual Information and Stochastic Modeling . 

The reflected energy from the target is detected by the scanner 
and corrupted by various noise sources. If S is the "noise-free" signal, 
Y the observation and N a random disturbance, then 

Y = S + N (48) 


the reduction of uncertainty about S obtained from Y is called the 
average or mutual information between the observation and original signal. 
Since the reconstruction of the reflected signal from the noisy observa- 
tion is the highly desirable capability, the comparison of such average 
information- and selection of these bands with the highest information 
content is chosen as- a means of spectral band selection. Let 


S = 
n 





s ) 
n 


and 


Y 

n 


(y- L » y 2 > 


*n> 


where s. and y. are the coefficients of the orthonormal (K-L) expansion 
i i 

of Y and S, then the mutual information between Y and S is given by [3]. 


I(Y,S) 


hlog 


det C 

. 1 

det C 


(49) 


where C and C are the covariance matrices of (y., i = 1, 2, ...) and 
y n i 

(n. “ No/2, i = 1, 2, ,..) and No/2 is the two sided spectral density of 
the additive white noise. Equivalently I(Y,S) can be represented in 
terms of the Wiener-Hopf optimum filter impulse response, 


I(Y,S) 



h(A, A) dA 


( 50 ) 



50 . 


h(A,A) provides an estimate of S from Y with a minimum mean-square error. 
This relationship, however, is not- a practical-method of evaluating I(Y,S) 
since the actual solution of the Wiener-Hopf integral itself is a nontrival 
task. This problem can be circumvented by a d-is.cret-e state variable 
formulation, i.e.,- 


s.Ok+1) = j>s(k) + FW(k) ke[A 1 ,X 2 ] (51) 


where 


s (k+1) 


Sl (k-fl) 

s 2 (k+l) 


s (k+1) 
n 


f is an (nxn) .matrix 

JT is an (nxl) vector 

W(k) = a descrete independent Ga.ussian 
zero mean random process 

The formulation of the problem in the discrete domain -provides a practical 
way of computing h(.A, X) through Kalman filtering techniques. The discrete 
version of Eq. (50) is given by 

I,(Y,S) * hi Mk„k) 

k e ("52) 

The discrete nature of this approach makes the evaluation of Eq. (52) 
considerably more practical than its continuous counterpart. This is due 
to the fact that the Wiener-Hopf equation is easily solved in only those 
cases for which the analytical form of K s (X,u).„ the signal covariance 
function, is fairly simple,, not likely for most -random processes encountered 
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in remote sensing. Since h(k,k) is dependent on the parameters of Eq. (51) 
a concise representation of _s(k+l) is needed. 

The general form of Eq. (51) is given as an autoregressive (AR) 

model 

m ! m 2' 

1 a.s(k-j) + J b.t|»(k-j> + W(k) (53) 

j=l J j=l J 

The spectral response at the discrete 
wavelength k. It is a Gaussian random 
process. 

zero mean independent Gaussian disturbance 
with variance p . 

deterministic trend term used to account 
for certain characteristics of the 
empirical data. 

are unknown constant coefficients to be 
determined. 

The order of the AR model. 

The identification selection and validation of general AR models 
for the representation of a random process is a well developed technique 
[l4,15]. The identification of an appropriate model provides the 
necessary parameters required for the evaluation of I(Y,S) in Eq. (52). 

The model selection process for a selected number of ground covers has 
been carried out [3] leading to the ranking of a set of spectral bands 
according to the criterion outlined previously. A summary of the 
experimental results are given below. 

Data Base and Model Selection . 

Two different sets of empirical data are used to demonstrate the 
techniques developed here. The first set consists of observations o,f 
wheat scenes. The second set consists of several vegetation cover 
types such as oats, barley, grass, etc. For each scene the spectral 
responses, collected by the Exotech 20C field spectroradiometer , are 
averaged over the ensemble. It is thought the resultant average 


js(k) = 
s(k) = 

w(k) = 
’f' (k- j ) = 
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;sp.ectral .response provides a relatively good -data set for demonstration 
-purposes. 'Figures 21 ,and 22 show -the average response for both cover 
types . 

in order that the study be carried out in a context that is rela- 
tively realistic for multispectral scanners,, the spectral response -of 
the two data -sets is -divided into a number -of spectral bands, The' division 
is fairly arbitrary, but each band must contain a sufficient number of 
-data points to ensure accurate parameter estimation for model identifi- 
cation. The spectrum is divided into 9 bands from 0.45 pm to 2.38 ym 
with two gaps in the 1.34-1.45 yim and 1.82-1.96 .ym range -due to atmos- 
pheric absorption, see Tables 5 and '6. 


Table 5.’ Spectral Bands for Wheat Scene. 


Band Number 

1 

2 

3 

4 

5 

6 

7 

8 
9 


Spectral Wavelength 
Interval (ym) 

0.-45-0.54 

.0.54-0.62 

0.62-0.71 

0.71-0.85 

0,. 85-0. 99 

0.99-0.13 

1.13-1.34 

1.45-1.82 

1.96-2.38 


Table 6. Spectral 
Band Number 

1 

2 

3 

4 

5 
■6 

7 

8 


■Bands for Combined Scene 

Spectral Wavelength 
•Interval (ym) 

0.45-0.54 

0.54-0. -62 

0.62-0.71 

0.71-0. .85 

•0 ..85-0.-98 

0.98-1.12 

1.12-1 ..30 

1.45-1.82 

1.96-2.38 


9 
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The next step is the identification and validation of models that 
would adequately describe the aforementioned spectral responses.. Three 
different models were tested and compared, (a) autoregressive (AR) , 

(b) autoregressive plus a constant trend (AR+C) , (c) integrated 
autoregressive (TAR) . following the standard selection and validation 
techniques, one of the above 3 -models is selected which describes the 
scene most satisfactorily. Tables? and 8 show the selected models for 
the wheat and combined scene respectively, IAR— 2 in Table 8 in a second, 
order IAR. 



Table 

7. Modeling of 

the Wheat Scene. 

Band 


Order of Model 

Type of Model 

1 


7 

AR 

2 


2 

AR 

3 


11 

IAR 

4 


1 

AR+C 

5 


1 

AR 

6 


2 

AR+C. 

7 


5 

IAR ' 

8 


8 

IAR 

9 


6 

IAR 


Table 8 

Modeling of the Combined' Scene . 

Band 


Order of Model 

Type of Model 

1 


11 

IAR- 2 

2 


2 

AR' 

3 


11 

IAR 

4 


1 

AR+C 

5 


3 

-AR 

6 


1 

AR 

7 


•9 

AR+C 

■8 


'8 

TAR 

9 


1 

AR 
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Spectral Band Selection . 

It was initially stated that the information content of a set of 

spectral bands can be used in the selection of an optimum subset. Here, 

the preceding regression analysis will be used to evaluate the mutual 
information between the reflected energy and the observed signal in the 
9 spectral bands under study. 

The first step is the computation of the average information in y(k) , 

the received spectral process, about s(k). The reflected spectral scene 

response as a function of spectral bandwidth for each band of both 

scene types. The average information is computed for several values of 

2 

the noise variance, a Appropriate software is developed to carry out 

the calculation of I(Y,S). Figures23 and 24 show .the variation of 

I(Y,S) in nats for the wheat and combined scene in band 1. Similar plots 

2 

are shown for the infrared band 7, Figures 25 and 26. Selecting a =10 
for demonstration purposes, the average information for wheat and combined 
scenes are tabulated in Tables9 and 10 for 9 spectral bands. 


Table 9. Average Information for Wheat Scene I 


Band 

1 

2 

3 

4 

5 

6 

7 

8 
9 


I (Y , S ) nats 
34.50 
10.52 
20.35 
30.00 
44.96 
37.20 
60.31 
34.80 
50.10 
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Table 10. Average Information for Combined Scene Band. 


Band 

1 

2 

3 

4 

5 

6 

7 

8 
9 


I'(Y,S) nats 
41.33 
16.17 
22.98 
40.08 
45.73 
40.96 
78.25 
64.15 
74.19 


Using the information content of each band as an optimality criterion, 
the 9 spectral bands can be ranked, see Table 11. 


Table 11. Order of Preference of Spectral Bands for the 
Wheat and Combined Scenes. 

Rank Wheat Scene Band Combined Scene Band 


1 

2 

3 

4 

5 

6 

7 

8 

9 


7 
9 

5 

6 

8 
1 
4 
3 
2 


7 
9 

8 

5 
1 

6 
4 
3 
2 


The top 6 highest ranked bands, although ordered differently, are the 
same for both cover types. Moreover, other than band 1 which is in the 
visible portion, the remaining 5 are all in the infrared portion of the 
spectrum. Thus, relative to the average information criterion, the 
infrared portions of the spectrum is generally preferred to the visible 
portion since bands 2 and 3 are ranked lowest for both the wheat scene 
and combined scene. 
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The selection of a subset of the available spectral bands using the 
idea of their information content is a new approach in band selection 
and requires further investigation to evaluate its optimality in more 
concrete terms. One' of- the most useful optimality criterion is the 
selection of these bands that maximize the overall classification accuracy. 
No -documented relationship exists between the average information contents 
of a set of bands and, the subsequent class separability. It is true 
however, that such information measure is directly related to the optimum 
Weiner filter thereby providing a basis for the optimality of this 
ranking technique. 

3. THE UNIFIED SCANNER ANALYSIS PACKAGE BLOCK DIAGRAM 

The identification and development of a set of individual techniques 
and algorithms is only the first step toward a complete system simulation 
package. The usefulness of this package is fully realized only when the 
elementary modules are interconnected in a logical and clear fashion. The 
objective here is the integration of the available processors such that 
starting with a raw data base, the question of optimum spectral bands, 

IFOV size and the noise model can be answered with the classification 
accuracy as a primary performance index. 

3.1 System Structure 

One realization of such simulation model was shown in Figure 4 and 
is repeated here for convenience in Figure 27. USAP is basically composed 
of three distinct parts (a) a spatial path, (b) a spectral path and (c) a 
means for classification performance estimation. In the following 
individual software modules are discussed. 

Classification Accuracy Estimators . 

There are two classification performance estimators available 
(a) the analytic classification accuracy predictors and (b) the stratified 
posterior performance estimator. 



Figure 27 


The Block Diagram of the Unified Scanner Analysis Package (USAP). 


LARSYS 


CORELAT 


SCANSTAT 



Specification 


o-> 

o 


Implies path for spectral/spatial statistics 
Implies path for data or spectral functions 


»x)> 
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Analytic Classification Accuracy Predictor . The ACAP algorithm 
discussed in Section 2.1 is the primary processor in evaluating the 
performance of a scanner system when the probability of correct classifi- 
cation is defined as the primary performance index. This piece of 
software, as shown in its theoretical development, requires one major 
input in the form of the population statistics. In order to facilitate 
the operations, the format of the statistics deck is chosen to be 
identical to the one produced by LARSYS statistics processor although 
it contains a considerable amount of redundancy such as field coordinates. 
These cards are skipped. Among other user— supplied information is the 
desired spectral bands to be used in the analysis and the ■ sampling and 
grid fineness in the form of number of elementary cells per axis. 

There is obviously a trade-off between the estimator's variance, a 
decreasing function of the grid size, and the computation time. If 
N is the number of spectral bands and n the number of cells per axis, 
the per class number of cells to be tested in a set of M quadratic 
discriminant functions is n^. This exponential relationship calls for 
a careful selection of n particularly for a high dimensional space. 

On the other hand a small variance is very much required property of any 
estimator. 

The relationship between grid structure and the estimator's variance 
has been covered in detail £2]* It was shown that the classification 
accuracy obtained using ACAP exhibits a relative independence from n for 
n > 12. This property is preceded by a fairly rapid rise to a steady 

*■** a 

state value after which the magnitude of the P^ oscillations is within 
0.5% of the true value or the Monte Carlo derived reference. The choice 
of n is ultimately decided by the user depending on his specific needs 
and after some experimentation. Initially, however, a default option of 
n = 9 cells per axis is considered to provide quick turn-around time 
while keeping the quality of the estimate high. The output, in addition 
to the classification accuracy estimate, contains information on the 
transformed class statistics, .feature space and sampling grid structure. 
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Stratified Posterior Classification Performance Estimator . This, 
is the software implementation of the' iaigofithm, discussed in Section 2.2. 
The maximum conditional aposteriori class, probability is the criterion 
for classification and error estimation purposes.. The program does, 
not provide any options and the size of the internally generated random 
data is fixed. ACAP and SPEST produce different, bpt very clos.e results. 

Spatial Path . 

Data Base. The input data to the spatial scanner model is via the 

4 ’ " 

multispectral image storage tape containing satellite or aircraft 
collected data. This tape has been reformatted and is compatible with 
any LARSYS processor. 

Data Retrieval. The individual software units can access the 
> * 

available data base through various system support routines or any of 
the LARSYS processors. 

Spatial Correlation Analyzer . The determination of the scanner 

characteristic function requires a knowledge of the spatial properties 

of the input data, therefore a class conditional estimate of the spatial 

auto and crosscorrelation functions is needed. Let f^(x,y) be a two 

dimensional image of size N x N pixels in the kth spectral band then 

o o 

the spatial autocorrelation function estimate is given by [16] . 


N -r N -n 
o o 


R kk (T?r,) = c J J £ f k (i ’ j) ” t f k (i + T ’ 3 + n) " 


i=l j=l 


x,n = 0, 1, . . , , n Q -l (54) 

where p, = E{f k (x,y)}. The multiplicative factor C can be chosen to. be 
one of the following 


C 


I 


1 

(N — r) (N -n) 
o o 


(55) 
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C 


2 



if is known and C = then E{R^^(T,n)} - h k t ^ ien 

neither selection of C or C 0 will produce an unbiased estimate. The 
actual derivation of the mean and variance of when the mean is 
estimated is rather complicated. The bias of the estimate in one 
dimension is given by [16], 


I , N - j x 1 . _n 

E{R, (t) - R(r) } = ^^-RCt) ° N Var {p k > + o(N q ) (56) 

o o 

From Eq. (56) it follows that R^x) is asymptotically unbiased. This 
result can be extended to the two dimensional functions provided the 
autocorrelation function is separable along each spatial axis . In 
general the maximum lag, n^, must be chosen such that n^ << N^. As a 
rule of thumb it is desirable to keep the maximum lag less than one- 
tenth the sample size N . This will tend to avoid certain instabilities 
that can occur in autocorrelation function estimates. The across-band 
correlation function estimate is obtained using an identical relationship 
to Eq. (54). 

The empirically obtained spatial correlation matrix needs further 
processing to be used in the scanner spatial model developed in Section 
2.3. Specifically a Markov correlation model is fitted to the experimen- 
tally obtained R kk (r,n). By invoking the separability assumption for 
small lag values, 


/V A A 

R(xyn) - R(t) R(n) (57) 

where no subscript indicates either auto or crosscorrelation function. 
Table 12 shows the magnitude of the errors involved in carrying out 
this approximation on an aircraft data set. The error is expressed as 
a percentage of the experimental values. An exponentially dropping 
function is then fitted to the individual correlation functions along 
the sample and line directions. The fitting Is accomplished using a 
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Table 12. Error Matrix for Cross Correlation Function Approximation 
Between Channels 2 and 8. 


1.00 

.92 

.81 

.69 

.59 

.50 

.44 

.93 

.88 

.78 

.67 

.56 

.48 

.41 

.73 

.71 

.64 

.54 

.44 

.36 

.3 

.48 

.47 

.43 

.36 

.28 

.21 

.16 

.30 

.31 

.29 

.24 

.18 

.12 

.08 

.23 

.25 

.24 

.21 

.16 

.12 

.08 

.22 

.24 

.24 

.22 

.19 

.15 

.12 


1.00 

.92 

.81 

.69 

.59 

.50 

.44 

.93 

.86 

«75 

.64 

.55 

.46 

.4 

.73 

.67 

.6 

.5 

.43 

.36 

.32 

.48 

.44 

.38 

.33 

.28 

.24. 

.71 

.30 

.27 

.24 

.2 

.17 

.15 

.13 

.23 

.21 

.18 

.16 

.13 

.11 

.1 

.22 

.2 

.18 

.15 

.13 

.11 

.1 

) 

0 

0 

0 

0 

0 

0 

) 

2.2 

3.8 

4.5 

1.8 

4.16 

2.4' 

) 

5.6* 

6.2 

7.4 

2.3 

0* 

6.7! 

> 

6.4 

11.6 

8.3 

0 

17.5 

23.8 

) 

13.0 

17.2 

16.6 

5.5 

20 

38.4 

> 

16 

25 

23.8 

18.7 

8.3 

20 

> 

16.6 

25 

31.8 

31.6 

26.6 

16.6 
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weighted least square approach to the logarithm of R(x) and R(r|) . The 
slope of this linear fit determines the adjacent sample or line 
correlation coefficients. Specifically let 


F(i) = In R(i) i = 0, 1, h q -1 


then the parameters of the linear fit, are given by [17 ] 


T -IT 
C = (H W H) H W F 


(58) 


where 


0, 1, . . . , n -1 
o 

n , n +1 , . . . , 2n -1 
o o ° 

and the diagonal weighting matrix, W 


H(i) =1 i = 

H(i) ■ 2n -1-i i = 
o 


(59) 


(n -1-i) 

W(i) = a ° i = 0, 1, ...» n Q -l (60) 

with a as the weighting matrix diagonal base, o < a < 1. The weighting 
matrix via the control parameter assigns a progressively smaller weight 
to R(x,ri) for succeeding lag values. This weighting is necessary since 
the properties of the correlation functions show an increasing deviation 
from the underlying assumptions of separability and Markovian structure 
for higher lag values . 

A complete specification of the spatial properties of the available 
spectral classes requires determination of 


N + 2! (N-2) ! 

auto and crosscorrelation functions per spatial axis per class where N 
is the total number of spectral bands used in the analysis. The imple- 
menting software contains various default provisions in case the corre- 
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lation properties of the input data differs considerably from the afore- 
mentioned assumptions. 

The user specifies the area to be correlated by a run table entry 

run number followed by the field coordinates. In order to perform either 

the auto or crosscorrelation operations, appropriate spectral bands (s) 

need to be specified. The maximum lag, n^, in computing R(r,ri) is also 

a variable and is entered as a percentage of the image size in pixels. 

The value of n is dependent on the size of the area to be correlated, 
o ^ 

Since the magnitude of R(r,ri) for Landsat data is generally negligible 
for more than 4 or 5 pixels lag, n^ as a percentage can take on small 
values for large fields and vice versa. 

'Scanner IFOV Model . The scanner IFOV software is the computer 
implementation of the scanner characteristic function discussed in 
Section 2.3. The input consists of (a) spectral covariance matrix, 

(b) spatial correlation matrix along the samples (c) spatial correlation 
matrix along the lines and (d) IFOV size in terms of the number of high 
resolution pixels. A standard LARSYS statistics deck produced by the 
statistics processor constitutes the first item. The spatial correla- 
tion information in entered through an K x N symmetric matrix the (i,j) 
element of which 



is the pixel-to-pixel correlation for bands i and j , a„ is estimated 
by the spatial correlation software. The IFOV size is expressed on a 
relative scale in terms of the number of high resolution pixels within 
1 IFOV of the scanner PSF, e.g., 1, 2, 3, etc. There are two choices 
available for the functional form of the PSF, Gaussian and rectangular. 


The output generated by this software module is a spectral statistics 
deck which is the input class statistics transformed by the scanner 
characteristic function. This deck is used as input to the ACAP processor. 
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Additive Noise Model . By virtue of the parametric approach adopted 
here, the incorporation of the noise effect takes on a simple form. The 
noise statistics is characterized by a zero mean vector and a diagonal 
covariance matrix with zero off diagonal elements- This matrix is then 
simply added to each class covariance matrix, 


1 =1 + 
^ ^ 

f f 


n„ 


n„ 




(62) 


The diagonal elements of £ , a , determines the SNR in each spectral 

— N n i 

band. By the appropriate selection of a 's, different SNR can be 

2 n - 

specified for each band. Let a., denote frhe variance of the noise-free 

f k . . 

signal at the scanner output, then the SNR in the kth spectral band is 
defined by 


(SNR) k = a\ io 2 n (63) 

k k 

The choice of equal of unequal SNR in different bands based on experimental 
or theoretical results is at the analyst's discretion. 

Spectral Path . 

Data Base and Retrieval . EXOSYS is a software package which provides 
access to field measurement data taken with a variety of field instruments. 
A brief overview of the EXOSYS package will be given here, with more detail 
available in the EXOSYS manual [18], Data is collected and stored, on 
magnetic tape in field measurements format. During the reformatting 
procedure the data is calibrated and ancillary information such as 
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weather readings soil conditions, and plant growth status is placed in 
the identification record for each run. 

There are three processors in EXOSYS which are used to access the 
field data information - IDLIST, GSPEG, and' DSEL. The IDLI-ST processor 
scans the tape and lists information from the identification record as 
required. One can use this information to select appropriate runs- to 
represent informational classes. 

The GSPEC processor provides a punched deck consisting of the 
numerical values of the spectral response function for all of the desired 
runs. One can select a set of run numbers as input and the output will 
consist of a punched deck. Options exist for plotting the spectral 
response functions for the desired runs. 

l 

The DSEL processor simulates rectangular spectral channels and 
uses data from the tape to evaluate the response in each channel for 
the ensemble. The inputs are the spectral band locations and the run 
numbers on the data tape.. The output is a set of statistics for the 
specified channels. 

Optimum Spectral Basis Function System . For the optimum spectral 
function calculation the output of the GSEC processor is required. The 
appropriate ensemble can be selected by specifying a set of identifica- 
tion parameters such as date of collection, scene type, run number, etc. 

The cards containing the numerical values for the spectral response 
functions are used and stored on- disk in a format which is compact and 
convenient for future processing by the program SPRDCT. The files that 
are stored on disk may be transferred to magnetic tape for future use to 
avoid repeating the procedure involving the EXOSYS package. SPRDCT 
requires some information to be entered at the terminal to provide ID 
information for the ensemble. A list of all runs used by run number is 
printed after the data has been stored on disk. 
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Optimum Basis Function Calculation . The calculation of the optimum 
set of basis functions for an ensemble is accomplished by solving the 
matrix equation 


<j>r = kw$ 


( 64 ) 


to get the eigenvalues a 2> . and the eigenvectors <f> , <f> 2 » <i> N . 

The matrix $ is the matrix of eigenvectors $ = <J> 2 , •••> <^1 and T is 

the diagonal matrix of eigenvalues 


a 


1 


r = 


o 


L 


0 


a 


2 



The matrix W is a diagonal matrix of weights 


W, 0 

0 w 2 


( 65 ) 


( 66 ) 


W 


N 


R is the covariance matrix for the ensemble. Let the mean vector for the 

tT 


ensemble be m = [m^, m^, m^] then 


R. . = E {(x.-m.) (x.-m.) } 
1 i 3 1 


( 67 ) 


The maximum likelihood estimate is 
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where N g is the number of sample functions in the ensemble. 

If we let A=KW, then A is a real general matrix. An algorithm for 
solving for the- eigenvalues- and eigenvectors of a real general matrix is 
available [19] and is used here, with only slight modification. The 
algorithm makes use of Householder's method and the QR double-step itera- 
tive process to compute the eigenvalues. The eigenvectors are obtained 
by the inverse iteration technique. A sorting routine was added to order 
the eigenvalues and the corresponding eigenvectors.. 

The required inputs are the data in> the appropriate format and the 
weight' function. The output of this processor is a set of N eigenvectors 
or bas'is functions punched onto cards. Also, the eigenvalues and means- 
square representation error are printed. The eigenvectors can.-be plotted 
using GCS subroutines. 

t 

Data Transformation and Statistics Calculation . The eigenvectors 
are used to perform a linear transformation on the original data vectors 
f X } . The transformed vectors {Y} have the desired -properties provided 
by the Karhunen-Lo&ve expansion. Each element of the transformed- vector 
is given by 


y± - ^ + *. 2 x 2 + . . . + 4, iN x n 


where <f>.- is the jth element of the ith eigenvector. 
3-3 


(69) 


The inputs to this processor ape the eigenvectors and the data- set 
stored on the disk. The output is the set of statistics for each class. 
The statistics are printed and punched on a deck of cards for future 
processing. 


4. USER'S OUIDE TO USAP 

i 

The block diagram of the scanner parameter study, Figure 27, is made 
operational by a collection of compatible software packages. Each module 
is individually compatible with the <LARSYS environment facilitating 
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incorporation of LARSYS processor in the overall system performance 
evaluation. This section provides a guide for the acquisition and 
execution of each program available on the LARS IBM 370/148 computer 
system. 

Prior to a discussion of the individual modules some general remarks 
are in. order. The access to this program library is simplified by the 
allocation of two special disk storage devices designated by DHSYS and 
DHDSK. The former contains the text version of the software while the 
latter holds the source. These devices can be accessed using the 
appropriate GETDISK commands. 

GETDISK DHSYS 


and 


GETDISH DHDSK 

these commands will establish the proper links in a Read-Only mode and 
USAP initialization is complete. In the following subsections the 
required input and necessary steps to run each program are discussed. 

4.1 The Classification Accuracy Estimators 

There arfe two parametric classification accuracy estimators 
available to the USAP user, (a) the analytic classification accuracy 
predictor (ACAP) and (b) the stratified posterior classification 
accuracy estimator (SPEST) . The theoretical aspects of these processors 
have been discussed in Section 2. Here is a guide to their software 
implementation. 

Analytic Classification- Accuracy Predictor 

This program evaluates the performance of a Bayes classifier when 
the populations statistics are multivariate normal. The following 
control cards are required. 
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.Control Word 'Description 

*ACAP This card specifies the particular processor 

requested. 


CHANNELS 


CLASSES 


GRID 


END 


The desired subset of the available channels 
is given here. Note that the numbers appearing 
on this card are the order of the selected 
channels not their actual number. For example, 
if the available channels are 8j 9, 12, 14 
and channels 8, 9 and 14 are requested CHANNEL 
card should read 1, 2, 4. 

This card specifies the name of each class. 

Each name must be placed in a field 7 charac- 
ters long followed by .a blank. The continua- 
•tion card, if required, must have the word 
"CLASSES" at the beginning. 

This quantity controls the quality of the 
estimate. The higher the number the closer 
the estimate is to the true Bayes error rate. 
(See 'Estimated CPU time" for more details.) 

This card signals the end of the control card. 
Stat deck follows immediately. 


Remarks . The ACAP processor in its current form is capable of 
handling up to 20 classes and 8 spectral bands. The extensions of 
these parameters presents no conceptual difficulty. The required statis- 
tics deck is a standard LAR.SYS produced deck with no modifications. It 
must be punched in the character format, however. 

How to Run the Program . Make sure the DHSYS disk has been accessed 
properly. One reader file consisting of the control cards followed by a 
statistics deck is required. Type ACAP in the CMS environment. Appro- 
priate terminal and printer output is produced. 

Example of control card set up 
*ACAP 

CHANNELS 1,3,4 

CLASSES SOYBEAN ALFALFA WHEAT 
END 


Since GRID card is not included, its default value (9)' is selected. 
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Estimated CPU Time . The execution time is quite sensitive to the 
GRID card specifying the number of cells per axis. For the default 
grid size and a 4-dimensional space it takes approximately 2 minutes of 
CPU time per class to provide the requested classification accuracy 
estimates. The CPU time is most sensitive to the dimensionality of the 
feature space. Hence if the number of spectral bands is limited (less 
than 4) considerable increase in GRID number is possible. The default 
number of cells per axis is considered to be the minimum while 'still 
providing acceptable performance. Increasing the parameter improves 
the quality of the estimate somewhat at the expense of higher CPU time. 
The choice is left at the user's discretion. 

Stratified Posterior Error Estimator. 


This program is identical in purpose but different in approach to 
the ACAP processor. Givaia set of multivariate normal populations, 

SPEST provides the classification accuracies associated with each class 
using an internally generated random data base. The different estimation 
procedures between the two methods is transparent to the user. 


Description of the Control Cards 


Control Word 
*SPEST 

CHANNELS 


CLASSES 


Description 

This card specifies the particular processor 
requested. 

The desired subset of the available channels 
is given here. Note that the numbers appearing 
on this card are the order of the selected 
channels not their actual number. 

This card specifies the name of each class. 

Each name must be placed in a field 7 charac- 
ters long followed by a blank. The continua- 
tion card, if required, must have the control 
word "CLASSES" in the beginning followed by 
the rest of the names. 


END 


End of the control cards. 
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Remarks . In usage, this program is identical to *ACAB. The 
standard. LARSYS statistics deck follows the control and disk immediately- 
Erinter output contains the estimated conditional classification, 
accuracies. By virtue of their separate, approaches, *ACAP and' *SPEST 
provide different, but very close, estimates of the correct -classif ication 

t 

accuracies. 

i 

How to Run the Program . The reader file contains the control cards. 

followed by the LARSYS statistics deck. A sample, control card deck 

follows: 

i 


*SPEST 

CHANNELS 1,.2,4 

CLASSES ALFALFA SOYBEAN WHEAT ' 
END 


4.2 Spatial Path 

The spatial path in USAP consists of two main; software units. The- 
spatial correlation analyzer, CORELAT and. the S.GANNER IF0V model,. 
SCANSTAT. 


Spatial Correlation Modeling . This .program is a. 2-dimensionali 
spatial correlator the primary output of which, is. a. normalized, spatial- 
auto (cross) correlation matrix for any specified; area- The user speci- 
fies the coordinates of the desired, segment in the form, of anj initial’ and 
final line and column along with, the appropriate spectral bands;(s) . 
Following the* estimation, of the. correlation- matrix,, th.e, exponential fit 
option, if invoked, will fit an exponentially dropping function, to- the 
experimental values' of or Ry^On)' using a weighted; linear least 

squares technique. 
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Description of the Control Cards 


Control Word 

Control Parameter 

Description 

* CORRELATE 


This card specifies the particular 
processor requested. 

INPUT 

RUN ( . ) 
TAPE ( . ) 
FILE(.) 

Run number of the desired area. 
Tape number of the desired area. 
File number of the desired area. 

BLOCK 

LINE ( . , . ) 
COLUMN ( . , . ) 

< - 

Initial and final lines.' 
Initial and final columns. 

FUNCTION 

AUTO 

CROSS 

Autocorrelation function requested. 
Crosscorrelation function requested. 

CHANNELS 


Channels used for correlation 
operation. 

SAMPLELAGt 


Maximum cross track lag used as 
a percentage of the total number 
of samples. 

LINELAGt 


Same as SAMPLELAG except for 
along track lag. 

EXP0FIT+ 


If included exponential fitting 
operation is carried out. 

END 


End of control cards. 


Remarks . This program is capable of processing areas containing 
up to 2400 pixels. The maximum lag default is set at 20 percent of 
the total number of lines and columns. Both quantities can be altered 
by user supplied control cards. The exponential fit option provides a 
pixel- to-pixel correlation coefficient for the channel(s) specified. 

This number is computed from the estimated parameters of the exponential 
correlation model. 

How to Run the Program . The only required reader file is the 
control card deck, an example of which follows : 


f optional 
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♦CORRELATE 

INPUT RUN (74028500), TAPE ('268.9), FILE (3) 
BLOCK LINE (1,25), COLUMN (1,25) 

FUNCTION AUTO 
CHANNELS 2 
SAMPLELAGt' 25 
LINELAGt 25 
EXPOFITf 
■END 


after DHSYS disk has been properly linked to, type CORELAT to start 

t 

execution. Appropriate terminal and printer output is .generated. 

Scanner IFOV Model . 

{ 

This program computes the spectral statistics of a population at 
the output of a multispectr’al scanner provided the data spatial correla- 
tion approximately follows a Markov model. The scanner IFOV shape is 
limited to either a Gaussian or rectangular shape. No assumptions 

a 

are made or indeed required about the type of the population statistical 
distribution. 


Description of the Control Cards 


Control Word 
* SCANS, TAT 

CHANNELS 


CLASSES 


IFOV 


Description 

This card specifies the .particular processor 
requested. 

The desired subset of the available channels 
is given here. Note that the numbers appearing 
on this card are the .order of the selected 
channels not their actual number. 

This card specifi.es the name of' each class. 

Each name must be placed in. a field 7 
characters long fo 1-1 owed by a blank. The 
continuation card, if required, must have 
the control work "CLASSES" in the beginning 
followed by the rest of the names.. 

This card specifies the -IFOV size of the 
scanner in .terms of -high resolution input 
pixels. 


t optional 
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APERTURE 


The choices here are "GAUSSIAN" or "RECTANGULAR." 


SNR+ 


Output signal energy to noise energy in dB. 


PUNCH+ The output statistics is punched out in an 

ACAP/SPEST format. Redundancies are added 
to replace field description cards. 

END End of control • cards . 


Remarks . This program is limited to 20 classes and 8 spectral 
bands. Execution time is quite short and extension of those parameters 
is straightforward. The input data immediately following the control 
cards consists of 3 parts: 

1. Standard LARSYS statistics ' deck in character format . 

2. Spatial correlation parameters (cross track) are entered via a 
NXN symmetric matrix where N is the number of channels . The 
(i,j) element of this matrix is the adjacent sample correlation 
between channels i and j . The lower triangular part of this 
matrix is punched in a 5 (E13.7,1X) format. 

3. Spatial correlation parameter matrix except for along track 
pixels. 

The above decks follow the control cards in the order listed. The 
signal-to-noise ratio is defined as the ratio of the output signal 
energy in a particular .channel (diagonal element of the class covariance 
matrix) to the noise energy in the same bands expressed in dB and defined 
by 


(SNR) - 10 log a 2 /a 1 

k k 


where k refers to the particular spectral band. 


t optional 
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How to Run the Program . One reader file consisting, of 4 consecutive 
decks and appropriate link to DHSYS disk is required before the program 
execution. An example of a control card, set up follows: 


*S CANS TAT 
CHANNELS 1,2,4 

CLASSES ALFALFA SOYBEAN WHEAT 
IFOV 2 

APERTURE GAUSSIAN 

SNRt 10 

PUNCH! 

ENDt 


The output consists of an ACAP/SPEST compatible statistics deck for 
the modified population. This deck can be used in the ACAP/'SPEST processors 
to obtain the new set of classification accuracy estimates. 

4 . 3 Spectral Path 


Tjlhe spectral path in USAP consists of three main pieces of software 
(a) data retrieval through EXOSYS processor (b) optimum spectral function 
calculation and (c) data transformation and statistics calculation. 

Procedure for Computing and Evaluating Optimum Basis Functions . 

Data Retrieval . The data retrieval system is stored on EXOSYSDV and 
it is necessary to define storage as 768K. A card file containing the 
data points for each run will be constructed on a temporary disk 
(25 cylinders) . It is desirable to make the temporary disk a P-disk 
and the permanent user disk a T-disk. 

In CMS 


1 nnt-lnnal 


RELEASE 191 P 
LOGIN 191 T 

GETDISK TEMP MODE P25CYL NO-UFD 



In CP 
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I EXOSYSDV 
CCINPUT TERMINAL 
RUN EXOSYSDV 


The control cards will be entered through the terminal. The tape on 
which the data is stored must be specified as well as the cover type 
and the collection date. A typical sequence of cards is as follows: 


$ TAPE 4896 
$ GSPEC 

GRAPH SPEC (SPRING WHEAT) , DACO (770508) 

LIST NO LIST 

OPTIONS PUNCH, NOGRAPH 

END 

$ END 
$ EXIT 


The runs taken over Williams Co., North Dakota for May 8, 1977 are on 
tape 4896. The crop species Is spring wheat and the collection date is 
May 8, 1977. The output is a deck of cards with one hundred data values 
for each run punched onto the cards. Header information must be flushed 
at the line printer. Return the system to CMS and read the cards onto 
a disk file SPR100 DATA. The number of records in the file is equal to 
the number of runs. This number should be recorded as it will be needed 
later. This procedure is repeated for the second class. The cards from 
GSPEC are read onto the file INPUT DATA and the number of records 
recorded. The two files are combined by typing (in CMS) 

COMBINE SPR100 DATA PI SPR 100 DATA PI INPUT DATA PI 

ERASE INPUT DATA 

This procedure is repeated until all classes have been included in the 
file SPR 100 DATA. The crop species SUMMER FALLOW and PASTURE are used 
in GRAPH SPEC(.) to complete the data set. 



At this point the program SPRDCT is loaded and run,. -;A disk file 
will be created using DSRN of 2 and file type ;EUNC. The following 
information is requested at the terminal by SPRDCT 


Experiment Number 


100 

Number of Classes 


3 

Number of Sample Points 

per Run 

•100 

(Dimensions) 



Class Name Wheat 

Fallow 

Pasture 

-Number of Samples 

per Class 664 

437 

16.4 


The information is requested and is entered -between the .slash 
marks, right justified. 

At this point the data is ready to be used by the .system. It is a 
good idea to store the file on tape for future use. Type 

TAPMOUNT 156 TAP 2 RI 
T DUMP SPR 100 FUNC PI 

The tape on which the data is stored is 156 in this example. 'To recall 
the data to the disk type 

TAPMOUNT 156 TAP2 R0 (I-f not already '.mounted-) 

T SLOAD SPR 100 EUNC 

Note that the P-disk must .be a large fairly empty -disk (10 cyl) 


The format for data storage is 
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ID 


Class 1 


Run 1 

100 Data 
Values 


Class 2 


Run 2 


100 Data 
Values 




ID Information Record 


Item Words 

Date 1-15 

Exp. Number 16 

Number of Classes 17 

Number of Dimensions 18 

Number of Samples for Class 1 21 

Number of Samples for Class 2 22 

Number of Samples for Class 3 23 

Number of Samples for Class 4 24 

Number of Samples for Class 5 25 

Number of Samples for Class 6 26 

Number of Samples for Class . 7 27 

Name of Class 1 30-39 

Name of Class 2 40-49 

Name of Class 3 50-59 

Name of Class 4 60-69 

Name* rvF HI ass 5 70—79 



Name of Class 6 
Name of Class 7 


80-89 

90-99 


Optimum Spectral Functions Calculations . 

Once the data set is on disk it is necessary to issue the following 
CMS commands to compute the eigenvectors. 

FILEDEF 2 DSK-P4 SPR100 FUNC RECFM VS LRECL 400 BLKS-IZE 
400 (PERM 

CP SP PUN TO USERID 
LOAD SPOPTM (XEQ 

4 . 4 Example Outputs 

This section presents a sample output for the individual software 
units used in USAP. The example set consists of sample outputs from the 
ACAP and SPEST processors. CORELAT and SCANSTAT in the spatial path and 
SPOPT and SPTES in the spectral path. Graham Co., Kansas is used as the 
test site. 


Classification Accuracy Estimators . 


The following control card set up is used for the ACAP processor and 
output is shown in Table 13. 


*ACAP 

CHANNELS 1,2, 3,4 

CLASSES BARESOI CORN SOYBEAN WHEAT 
END 


The required control cards for the SPEST processor are 
*SPEST 

CHANNELS 1, 2,3,4 

CLASSES BARESOI CORN SOYBEAN WHEAT 
END 


the output is shown in Table 14. 



Table 13. *ACAP Sample Output 


ANALYTIC CLASSIFICATION ACCURACY PREDICTION - A C A P 


CLASS WHEAT 


SAMPLING GRID CHARACTER! ST I CSO 
GRID SIZE” 9 CELLS PER DIMENSION 
TOTAL NO OF CELLS IN THE GR 10” 6561 


TRANSFORMED FEATURE SPACE CHARACTERISTICS!) 
EIGENVALUESO 

2.9182E 01 1.14I3E 01 1.3430E 00 

EIGENVEC TORSO 


2. 1045E-OI 
4.5844E-01 
-6.7334E-03 
8.63426-01 


2.0478E-01 

8.2659E-01 

1.9331E-01 

-4.8728E-0L 


TRANSFORMED MEAN VECTDRSO 


1.1654E 01 
-2.2766E 00 
1.3111E 01 
0.0 


2.3928E 01 
1.5Z4BE 01 
5.1333E 00 
0.0 


8.2319E-01 

-1.52396-01 

-5.3271E-01 

-1.2389E-01 


3.82006-01 

1.3871E-01 

-5.9814E-0I 

0.0 


7.61256-01 


4.85936-01 
-2.8873E— 01 
8.2390E-01 
4. 12836-02 


-5.4750E-02 
1.0838E-01 
1.0702E 00 
0.0 


Q o 

tv 

ll 


<© 




c? 

tr< £> 

^ & 


PROBABILITY OF CORRECT CLASSIFICATION* 93.486 I 


00 

u> 


***** total PR08 OF CORRECT CLASSIF ICAT IQH*****= 89.287 PERCENT 



Table 14. *SPEST Sample Output 


STRATIFIED POSTERIOR error estimator 
CLASS BARESOI 

probability of correct clas sif 

CC ASS CORN' 

PROB, ABILI. TV OF CORRECT CLAS SIF 

CLASS SOYBEAN 

PROSABI'LITY. OF CORRECT CLAS SIF 

CLASS WHE.AT' 

PROBABILITY; OF CORRECT C L A’ S- S I F 

OVERALL PROBABILITY: OF CORRECT RECOGNITION - 90;228 


I C A T I 0 N - 79.2181 


C A' T I 0 N » 92.7211 


CATION* 96.3591 


C A T I ON* 92.61'M 


<x. 
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Spatial Path. 


The control card set up for CORELAT is as follows: 


*CORRELAT 

INPUT RUN (74028500), TAPE(2689), EILE(3) 
BLOCK LINE(50, 98) , C0LUMN(50,98) 

FUNCTION AUTO 
CHANNELS 1 
EXPOFIT 
END 


The sample output is shown in Table 15. 


The control card set up for SCANSTAT is as follows: 


*SCANSTAT 
CHANNELS 1, 2,3,4 

CLASSES BARES 01 CORN SOYBEAN WHEAT 
IFOV 2 

APERTURE GAUSSIAN 
END 


The sample output is shown in Table 16. 

Spectral Path . 

The optimum basis function calculation and computation of the 
transformed data statistics comprises the main spectral processors. The 
sequence of required commands has been shown in Section 4.3. The 
following example is a weighted basis function calculation. 

The weighting function w(X) is zero for the water absorption bands 
near 1.4 and 1.8 micrometer and zero elsewhere on the interval (.4-2.4 pm). 
The printer output is shown in Table 17, listing the first 30 eigenvalues. 
The first 4 eigenvectors are sent in card format to the reader. They 
can be stored on the disk by issuing the command. 


0 READ EIG100 DATA T1 



Table 15. *C0RRELAT Sample Output 


TWO DIMENSIONAL SPATIAL CORRELATION' ANALYSIS 
CHANNELSO l I 


2-0 SPATIAL CORRELATION MATRIX 


1 .00' 

0..75' 

0.50 

0 .,70 

0 .,60; 0.44 

0.50 

0.45. 

0.36 

0;.38 

0.37' 0 .,3.1 

0.31- 

0.32 

0.27 

0‘.2S 

0.26 

0.23 

Ot.20' 

0.20' 

O’.lB. 

0..14 

0..15 

0,'. 14. 

0.10 

0.1-1', 0.12 


0.36 

O'. 2 9 

O'. 23. 

0. 33 

0.26', 

O'. 20, 

0.27 

0.21, 

0 .17 

0,. 24' 

O'. 1.7, 

0.1-3 

0.21, 

0.14 

0.1,1- 

0..16 

0.11 

0.09 

0. 13‘ 

0..10- 

0.07-- 

0.10-, 

0.07 

0.C5; 

o,. n 

0.09.' 

o-. 06; 


0.19' 

0. 13' 

0.06 

0.17 

Ov 12 

0..07 

0_.,l,4 

0.10- 

0.05- 

0,. 11 

o. os; 

0 .03'. 

0. 10 

0.08, 

0.04, 

0.07- 

0.05 

0 .02’ 

o;.,06' 

0.05- 

0,.03t 

0.06' 

Oi 04. 

0 .02,- 

0.031 

O’. 02 

0..00- 


INFORMATION'.'-’- - 


WEIGHTED?, LEAST- SQUARES FIT’ 


WE IGHTING MATRIX. DIAGONAL. BASE=0.,40’ 

WE.IGHTED' LSF. ERROR. {CROSS TR A CK ) =-, 0. 103756.9E-02,’- 
AOJACENT SAMPL-E CORRELATE 0N= 0.7 1.9 37.1 6 E 00,. 

WE IGHTEDL-SF; ERRORU ALONG! TRACK-I = , 0.6734634Er05- 
ADJACENT; L.INE CORRELATION-. 0.7G26460E 00: 



Table 16. *SCANSTAT Sample Output 


SCANNER OUTPUT STATISTICS 

APERTUREO GAUSSIAN 

IFOV SIZEO Z HIGH RESOLUTION PIXELS 

CLASS CORN 


INPUT COVARIANCE MATRIX 


9.29 




12.26 

19.79 



10.63 

16.37 

16.09 


A. 43 

7.14 

6.24 

3.45 


OUTPUT 

COVARIANCE 

MATRIX 


4.86 




4.61 

10.37 



4.00 

6.15 

8.43 


1 .67 

2.66 

2.34 

1.81 


O Q 

^ fed 

TJ Q 


W 

I 


> 

a 

c?d 


CO 


00 

'-J 
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LARQPATOkY for applications of remote sensing 
' PUPf.UF_ UNIVERSITY • 


SAMPLE FUNCTION INFORMATION 10 JULY, 1978 

EXP. NO 100 

NUMBER OF CLASSES 3 

class: WHEAT 

NUMBER OF SAMPLE FUNCTIONS. 664 

CLASS, FALLOW 

NUMBER OF SAMPLE FUNCTION^ 437 

CLASS PASTURE 

NUMBER OF SAMPlE FUNCTIONS 164 


WEIGHTING FUNCTION NUMBER 3 


N 

EIGENVALUE 

VAR (GAM) 

1 

311.4133 

307.0752 

2 

229.3520 

166.5619 

3 

21.1702 

1.4191 

4 

15.4660 

0.7574 

5 

8.8838 

0.2499 

6 

5.7765 

0.1057 

7 

3.5611 

0.0402 

8 

2.6887 

0 .0229 

Q 

2.3128 

0.0169 

10 

1.8363 

0.0107 

1 1 

1.41 95 

0 .0064 

1 2 

1.2304 

0 . 00*6 

n 

0.9300 

0 .0027 

I A 

0.6806 

0.0015 

15 

0.5217 

0 .0009 

16 

0.3517 

0.0004 

17 

0.3138 

0 .0003 

13 

0.2945 

0.0003 

19 

0.2771 

0 .0002 

20 

0.2336 

0.0002 

21 

0.2102 

0.0001 

22 

0.2092 

0 - 0 0 0 1 

23 

0.1792 

0.0001 

?4 

0.1514 

0.0001 

25 

0.1452 

0 .0001 

26 

0.1137 

0.0000 

27 

0.0976 

0.0000 

2P 

0.0851 

0 .0000 

2 Q 

0.0744 

n . oooo 

30 

0.0583 

0.0000 


V-AR (PHI ) 

mean-square 

ERROR 

0.008b 

299.156749 


0.0087 

69.804748 


0.0104 

4 f . 6 34554 


0.0126 

33. 1-68344 


0.0106 

24.204784 


0.0129 

1 o . 3 0b-26>9 


0 . 0c j7 

1^.947187 


0. 0573 

lc.2 5 b 4 8*9- 


0.063b 

9.945677 


0.0^44 

o. 10 5-3-93 


0.0663 

p . 6d98b0 


0 . 0o43 

5 • 45*494 


0.0336 

4 .52*464 


0.0338 

3.8^33-23 


0.0351 

3. 327086 


0.1 2o6 

2 ./ 975436 


0 . 3h 14 

2 . 6 o 1 5 9 9 


0 • 47o2 

2.3o7111 


0.3430 

2.090019 


0.21*6 

1 . 8 D 64-1 1 


35 . 72b2 

1.646236 


35.7220 

1.4.37049 


0.. 14.69 

1.237606 


0.5298 

1.106359 


0 . 52o6 

0. 9611.20 


0.0943 

<J . -8 h 7-453 


0.1201 

0 . 7498 54 


0. 1318 

0 * 664785 

V 

0.1133 

0 . 5-9 0.3 91 

S' ' 

0.R213 

0.5321 1-2 




-v ^ 


*4 

e 


Table 17. Eigenvalue and Mean- Square Representation Error for the Data 
Set. 
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A GCS routine was used to plot the graphs of the eigenvectors the 
first of which are displayed in Figures 28 through 31. 

The statistics for the first 4 terms are computed by using the 
same FILEDEF command as above plus 

CP SP PUN TO USERID 
0 PUNCH EIG100 DATA TI 
LOAD SPTES (XEQ) 

The program will ask 'NUMBER OF TERMS?', to get all 4 terms type '4.' 

The output will be a statistics deck with the following format: 

Card 

1 Number of Classes, Number of Terms 

2 Apriori Probabilities for each Class = 1/Number of Classes 
Mean Class 1 [Format (20A4) 1 

Covariance Matrix Class 2 
Mean Class 2 


The covariance are in upper triangular form. This statistics deck 
can be used as the input to the classification error estimator algorithms. 
Table 18 is a sample of the statistics obtained from the data set using 
the first 4 optimum basis functions. Also, the statistics were used as 
input to the classification performance estimator *SPEST to get an 
estimate of the probability of correct classification. 







MEAN VECTOR 

-200 .3751 22c 0 529 4.5466 >18.8553 

COVARIANCE MATRIX 
312.5391 

-24.4102 62.8445 

8.4233 -20.7412 14.0594 

-5.2539 -8.6057 4.5237 12.26-1-2 


MEAN VECTOR 

-202.1029 35.2806 5.7945 16.3977 

COVARIANCE MATRIX 
244.7227 

-69.3594 152.5959 

-15.8333 10.3513 24.0876 

-2.2070 1.3667 -0.7130 12.3877 


MEAN VECTOR 

-187.5431 54.8315 8.2578 19.0705 

COVARIANCE MATRIX 
286.1719 

-1.2813 168.4688 

-19.2971 -47.7388 29.0520 

10.3206 54.3401 -16.4006 26.-1763 

PROBABILITY OF CORRECT CLASSIFICATION FOR CLASS 1 = 0.9187 

PROBABILITY OF CORRECT CLASSIFICATION FOR CLASS 2 = 0.6624 

PROBABILITY OF CORRECT CLASSIFICATION FOR CLASS 3 = 0.9003 

OVERALL PROBABILITY OF CORRECT RECOGNITION - 0.-8270 


Table 18. SPTES and SPEST Sample Output Using the Firs't 4 Eigenvectors 
and Estimates of the Classification Accuracy. 
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5. SUMMARY 

The task of evaluating the performance of a multispectral scanner 
system while incorporating every spatial, spectral, electronic and tele- 
metric parameter is exceedingly complicated. The primary objective of 
this project has been the investigation of the two important, and most 
relevant in remote sensing applications, of scanner parameters, namely 
spatial and spectral. The development of analytic techniques for system 
performance evaluation differentiates the approach adopted here from 
other experimental methods. This property provides an ease of parameter 
manipulation not available through some heavily data dependent algorithms. 
Although the development of individual components of such systems is 
fundamental to the overall system operation, it is the logical and proper 
integration of individual modules that determines its ultimate processing 
capabilities. 

The Unified Scanner Analysis Package (USAP) is a fully integrated 
system with complete input-output compatibility of software units. It 
consists of a spatial and spectral path plus a shared unit providing the 
desired performance index in the form of probabilities of correct classi- 
fication. 

5.1 Classification Error Estimators 

The primary performance index throughout this study is defined as 

the probability of correct classification of the various populations 

present in a data set. In keeping with the underlying requirement of a 

parametric approach, the available LARSYS and other classification 

accuracy estimators using a randomly generated data base were deemed 

less than satisfactory. It is well known that the exact probability of 

correct classification is a multiple integral over an appropriate domain. 

The direct evaluation of such integral in a continuous N dimensional 

» 

space is a complicated and mathematically cumbersome task. This problem 
is circumvented by a deterministic sampling algorithm of the feature 
space preceded by an orthogonal transformation. This transformation 
when applied to the Gaussian probability density function of a particular 
class under consideration would conditionally decouple the feature space 
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and hence reduce the N dimensional error integral to a product of N one 
dimensional integrals each of which is a widely tabulated quantity. This 
algorithm requires the population statistics as the only major input 
and provides classification error estimates of high quality without 
excessively fine feature space quantization. The second classification 
accuracy estimator uses the maximum a posteriori principle coupled with 
a Monte-Carlo type integration technique. Although this algorithm is in 
a way dependent on a simulated data base, from a userspoint of view the 
difference between ACAP and SPEST are essentially transparent since both 
methods require the spectral statistics of the populations as their 
primary input. The aforementioned classification accuracy estimation 
techniques provide the basic tools for the scanner system performance 
evaluation. 

5.2 Scanner Spatial Parameters Selection 

The scanner spatial modeling algorithm and software consists of one 
main plus two supporting routines, i.e., IFOV modeling (SCANSTAT) , spatial 
correlation analyzer (CORELAT) and classification error estimator (ACAP) . 
The objectives of an analytical representation of scanner TFOV model 
is the, establishment of a parametric relationship between the system's 
input -and output statistics in terms of the class conditional mean vectors 
and covariance matrices. This relationship is -established using linear 
system analysis techniques extended to a 2 dimensional space. In order 
to derive any specific results two basic characteristics need to be 
specified: (a) scanner PSF and (b) ground scene spatial correlation 

model. 

The choice of a Gaussian shaped P.SF has been widespread in the field 
of image processing as applied to the Landsat data. This model closely 
approximates the averaging property of the -scanner aperture. An added 
feature of a Gaussian shaped PSF is the simplification of an otherwise 
intractable and cumbersome mathematics. Generally speaking the amount 
of information available about the spatial correlation properties of 
remotely sensed data is sparse. It has been frequently observed however, 
that the ground scene spatial correlation model approximately follows an 
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exponentially dropping function [2]. On the basis of previous experimen- 
tal evidence and the mathematical simplicity afforded by these assumptions, 
a Gaussian PSF and a Markov scene spatial correlation model is adopted. 
Like many other instances, the choice of the problem assumptions does not 
necessarily rest on their strict validity but also on the tractiability 
of the ensuing algebra. It is entirely conceivable that much more 
elaborate scene correlation models and PSF shapes can be envisioned. 

This approach, however, could and would complicate the underlying mathe- 
matics to the point where the gains initially expected from the more 
accurate model are balanced out. For simulation purposes the entire 
analysis is repeated for a rectangular shaped PSF although no currently 
operational Landsat is equipped with such a scanner system. 

Based on the foregoing discussion the scanner characteristics 
function is derived in a closed form. This function relates the input 
and output statistics as a function of the IFOV size and pixel- to -pixel 
correlation. SCANSTAT is the software implementation of this linear 
transformation. The auxilliary program, CORELAT, estimates the class 
conditional correlation functions and provides the best exponential 
curve fit to the experimental data using a weighted least-squares fit 
algorithm. The resulting output statistics is modified by additional 
random noise the power of which is computed from the specified SNR. The 
ACAP or SPEST processors provide the new classification accuracy sets . 

The probabilities of correct classification at the scanner output provide 
the basic information needed to evaluate the system performance under 
various operating conditions. For test purposes a hypothetical set 
consisting of 3 populations is selected and their statistics (mean 
vectors and covariance matrices) specified. The scanner output statist- 
tics and associated classification accuracies are computed for various 
IFOV sizes and scene correlations. The results are in close agreement 
with the numerically oriented experiments. For any fixed scene correla- 
tion, the population separability and hence the overall classification 
accuracy increases monotonically with IFOV size. The rate of increase, 
however, is a function of the scene spatial correlation. The classifi- 
cation accuracy increase per IFOV step is small for a highly correlated 
scene compared to a scene with a small adjacent sample correlation. 
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This property steins from the features of the 'scanner characteristic 
function and its particular weighting process- The addition of white 
Gaussian noise predictably degrades the output separability. The 
experimental results show that for a fixed IFOV size, SNR and classifi- 
cation accuracy increase monotonically . Same relationship exists between 
the classification accuracy and IFOV size when SNR is fixed. 

5,3 Scanner Spectral Parameters Selection 

The task of information extraction from remotely sensed data here 
primarily deals with the development of methods and techniques to 
select a set of spectral bands to enhance population separability. The 
first criterion employed in selecting a set of optimum spectral channels 
is the Karhunen-Loeve expansion of the ensemble of spectral responses 
associated with a cover type. This expansion provides a "set of optimum 
basis functions, a linear combination of which reconstructs the original 
stochastic process with a minimum mean square error. These basis functions 
in effect define a set of optimum windows in the electromagnetic spectrum. 
The associated software consists of EXOSYS data retrieval package, SPOPT 
spectral function claculation and SPTES data transformation and statistics 
calculation. The classification accuracy estimates used to check the 
resulting separability is obtained using either the ACAP or SPEST 'proces- 
sors . 


The se’cond approach employes an information theoretic concept for 
the specification of the optimal spectral bands. The observed spectral 
random process is modeled as the sum of a noise free signal plus an 
additive random noise component. For a candidate set of channels the 
quantity of interest, mutual information between the reflected and 
observed energy, is computed. The method consists of representing 
each random process as an autoregressive model. This type of represen- 
tation facilitates the evaluation of the mutual information when expressed 
in terms of the Wiener-Hopf filter PSF. Experimental results consists 
of selecting a wheat scene and dividing the continuum of electro magnetic 
spectrum into 9 distinct bands. In each band a proper autoregressive 
model is fitted to the particular random process. Following the 



97 


estimation of the parameters of the regression models, the average 
information content in each band is computed and on this basis spectral 
channels are ranked. Therefore on the basis of maximum mutual information 
optimality criterion, the top N bands represent the N "best" choice. The 
significant result obtained from the ranking method is that out of the 
6 top ranked channels 5 lie in the infrared portion of the spectrum 
thus future scanner systems used in remote sensing application should 
contain more infrared spectral bands according to this analysis. 

5.4 Conclusions 


This report has presented a brief description of the algorithms and 
results of the Scanner Parameter Selection culminating in the development 
of a Unified Scanner Analysis Package by proper integration of the 
available software modules. Although this report is the final document 
in this project, it actually represents the first step toward a well 
coordinated scanner system parameter study technique. The current struc- 
ture of USAP basically represents a skeleton of the future analysis 
packages. There exists a considerable software and theory development 
potential. The software by and large can take most of the streamlining 
to further facilitate their usage. Specific topics include extended 
diagnostic handling and error recovery capability, accelerated algorithms 
to further reduce execution times, etc. 

An overall * evaluation of the methods and results presented in this 
report shows that the objectives initially outlined have been success- 
fully met. The resulting analysis package, starting from a data base, 
produces specific guidelines on the selection of spatial and spectral 
parameters of a multispectral scanner system and it does so on an entirely 
analytic basis. In closing it should be pointed out that USAP can have 
a pivotal role- in any follow up project providing by far the widest and 
most economical parameter manipulation scope, fully complimenting any 
numerical or experimental scanner analysis techniques. 
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FILE. . . ACAP FORTRAN B1 

********** a** ************************ ***♦ **************** ************* 
SUPERVISOR TOR THE ANALYTIC CLASSIFICATION ACCURACY PREDICTION 


WRITTEN 08/16/78 


B I JAN G. HOBASSERI 


******** ********************** ******************************* ********* 

DESCRIPTION AND PURPOSE 

THE ANALYTIC CLASSIFICATION ACCURACY PREDICTION - ACAP - IS A 
SOFTWARE PACKAGE, THF PRIMARY OUTPUT OF WHICH IS A SET OF 
PROBABILITIES (IF CORRECT CLASSIFICATION FOR VARIOUS NUMBER OF 
CLASSES IDENTIFIED IN THE DATA SET. 

THE ALGORITHM IMPLEMENTED HERE DIFFERS FROM THE DATA DEPENDENT 
TECHY I UUfcS AND RATIO ESTIMATORS SUCH AS *CL ASS I F YPOI NT OF EAR SYS. 
WHEREAS LARSYS RETURNS TO THE DATA BASE FOLLOWING THE ESTIMATION 
OF THE CLASS STATISTICS IN ORDER TO PROVIDE THE CLASSIFICATION 
ACCURACIES, ACAP BYPASSES THIS STEP AND PROVIDES THE SAME 
QUANTITY DIRECTLY FROM THE STATISTICS DECK GENERATED BY * ST AT. TH 
METHOD EMPLOYED IS A MATHEMATICAL ONE (AS OPPOSED TO STATISTICAL! 
AN J INVOLVES PERFORMING A MULTIDIMENSIONAL INTEGRATION OF EACH 
PAD 3A IL t FY DENSITY FUNCTION OVER THE HYPERVOLUME OF CORRECT 
DECISION DOMAIN. 

DESCRIPTION OF CONTROL CARDS 

♦ ACAP 

THIS CARO SPECIFIES THE PARTICULAR PROCESSOR REQUESTED 
CHANNELS 

THE DESIRED SUBSET CF THE AVAILABLE CHANNELS IS GIVEN HERE. . 

IT IS IMPORTANT TO REMEMBER THAT THE NUMBERS APPEARING CN THIS 

Card is the urder of the selected channels not their actual 

NUMBER. FOR EXAMPLE IF THE AVAJLBLE CHANNELS ARE 8,9,12,14 

AND CHANNELS 8,9 AND 1% ARE REQUESTED, SUBCHANNELS CARD SHOULD 

READ 1,2. A. 

CLASSES 

THIS CARD SPECIFIES THE NAME OF EACH CLASS. EACH NAME MUST BE 

PLACED IN A FIELD I CHARACTERS LONG FOLLOWED BY A BLANK. THE 

C ON T IUE A 1 1 ON CARO , I F REOIUREO , MUST HAVE THE WORD 'CLASSES* 

IN THE BEGINNING rOLLOWfO BY THE REST OF THF NAMES. 


THIS QUANTITY CONTROLS THE QUALITY OF THE ESTI MATE. THE HlCHEi 
THE NUMBER THC BETTER THE ESTIMATE (ON A MACRO SCALE!. THE CPi 
TIME .HUWEVFR if VERY SENSITIVE TO THIS QUANTITY AND FOR 
INITIAL USAGE if 15 STRONGLY RFC CHHENDEO TO LEAVE (HIS TO ITS 
DEFAULT VALF OF 9 r.CLLS/aXIS UNTIL THE USER HAS ACQUIREG A 
FEEL FOR THE PROCESS. 


THIS CARO SIGNALS THE END OF THE CONTROL CAROS. STAT DECK 
FOLLOWS IMMEDIATELY. 

**• STAT OECK MUST BE PUNCHED IN CHARACTER FORMAT *** 

REMARKS 

THE PROGRAM IS CURRENTLY CAPABLE OF HANDLING 20 CLASSES AMO , 
8 SPECTRAL BANDS. THE EXTENSION OF THESE PARAMETERS IS ONLY UNI 
BY COMPUTATION TIME. THE USE OF MORE THAN 5 FE ATUR6S I S NOT 
RECOMMENDED UNLESS THE GRID SUE IS REDUCEO. TRIAL AND ERROR IS 
PA 08 ABLY THE QUICKEST WAY OF FINDING THE PROPER TRADE OFF. 

HOW TO RUN THE PROGRAM 

THE FOLLOWING EXEC FILE .NAMED ACAP, MUST RESIDE ON THE USER • S 


FILE. . . 


ACAP FORTRAN fll 


PRIMARY DISK 
ETYPEOUT OFF 

GLOBAL TXILIB SYSLIB CHSUB SSPLIB 
GETDISK LARSYS 
CP C C 

CP SP RDR HOLD 

L3AD ACAP BCOVAL (CLEAR NOMAP XEQI 
£E XI T 


AFTER THC- READER FILE IS LOADED , TYPE ACAP IN CMS ENVIRONMENT. 
THIS COMMAND INITIATES THE EXECUTION FOLLOWED BY APPROPRIATE 
TERMINAL MESSAGES AND PRINTER OUTPUT. 

EXAMPLE OP THE CUNTRCL CARO SET UP 


C *AC AP 

C CHANNELS 1,3.4 

C CLASSES BARESOI CORN PASTURE WHEAT 

C GUO 9 

C END 

8 DEFAULT VALUE OF GRtD IS SELECTED IF GRID CARO IS NOT 

C INCLUDED, 

C ESTIMATED CPU TIME 
C 

C WITH THE DEFAULT GRID SIZE AND USING 4 FEATURES, IT TAKES 

C APPROXIMATELY 2 MINUTES CF CPU TIME PER CLASS TC PROVIDE THE 

C REQUESTED CLASSI F ICA TIDN ACCURACY ESTIMATES. 

C TIE CPU TIME IS MOST SENSITIVE TO THE DIMENSIONALITY OF 

C THE FEATRE SPACE. HENCE IF THE NO OF SPECTRAL BANDS IS 

C LIMITED (LESS THAN 41 CONS I LERARLF INCREASE IN GRID NUMBER 

C IS FEASIBLE. EXPERIENCE HAS SHOWN HOWEVER, THAT THERE IS ONLY 

C A MARGINAL IFRaCTIuN UF PERCENT! I MPPROVFMENT I A THE 

C CLASSIFICATION ACCURACY ESTIMATE BY INCREASING THE GRID 

C SIZE BEYOND ITS DEFAULT VALUE. 

C 

C******** ************************************************ a ** a *• a*** * o* *> 
C 

REALMS A I 720 1 ,AA|T20) .SIGMA I 12B0 I ,M&[ I6D ! . HAAf 160) , M U60 i ,SD( 8 1 , 

1 DELTAIB) ,Pt 160 >,PC ( 20 1 , V 1 20 I , WV1 1 0 1 . WV2 1 8 I , XPR IN I 64) , 

2 RIGS! ,0 164) ,PR 1 1 31 ,EV (81 ,LLI8! ,UL I 8 5, HI 20 l,Z< 104 I ,DETI20I 

3 • E M A ( 1 60 1 « fc A I 72 Gl 
RE4L*4 ICSEH9C) 

INTEGER* 4 NADR (720), IC(20 ! , I NDX 1 20 1 .HEAD ( 20) ,SQIW 
I NT FGER* 4 LI$T(5!/»*ALA», ‘CHAN*. *CLAS», 'GRID*, 'EM) '/,1VEC(II. 

1 IC ARO 1 20 > 

INI E5CS*4 BLANK/ 1 »/,F STC RD/ ' L ARS 1 / 

INI EGER* 2 ICSELOOI.NCIIO), NCC I 30 ) 
f LOGICAL*! FLAG!5!/5*. FALSE./ 

COMMON /CAPCOH/ ND.NO IH, N S . NCL S , NP, I S f ZE , KNT R 
c COMMON /ACUCOH/ NT S , N TM, N SS , NS A , NSH , SD t M , NS I ZE , I SKP 

II FORMAT!//! 

IOIO FDRMATt IX, 'ERROR INCTLWRO. EXCCUTION TERMINATED.'! 
toil FORMAT IX, 'ERROR IN IVAL. EXECUTION TERMINATED.*! 

1012 f UR MAfilX, 'ERROR IN CHANEL. EXECUTIO,! TERMINATED,') 

1013 FORHATUX, 'MISSING CONTROL WORD. EXECUTION TERM! MATEO.* I 

1014 FORMA Ti l X. 'ACAP PROCESSING STARTEO'I 
£015 FOR MA 1 1 1 X , ' ALL CONIROL CAROS HAVE BEEN READ' ! 


C******** ********* ************* ******************* ***************** 

c 

C DECODE CONTROL CAROS 

C 

C**«* ****** ******* ********* ******************* ******** ,***•»***«*♦* 

c 

C SET THE DEFAULT GRtD SIZE 


DO 777 1*1,30 
ICSEL 111*0 

ncT i]»o 
CONTINUE 


ORIGINAL PAGE IS 

OF POOR QllALISY 

• *■ 



ACAP 


FORTRAN 01 


DO 778 1*1,90 
ICSETUI*- 50000.0 

FLAG141* .TRUE. 

CONTINUE 

LSZ-5 

1ER*0 

1NRU-5 

CALL CTLHRO nCARD,ICOL.ltST,LSZ,ieODE,INRO,IER) 

IFIIER.NE.O) GO TO 1001 

GO TO (99,101(102,103(104}. ICOOE 

FLAG! II*. TRUE. 

GO TO 100 

CHANNELS CARD 

CALL CHANEL ( I CARD ( I C OL (NCR. ICSEL, ICSET ,NCC, 900 1 
FLAG! 2 1 * • TRUE . 

ND1H-NCA 
GO TO 100 

CLASS NAMES CARD 

00 10 1*1.20 
HEAD! t ) » ICARDC 1 1 
FL AG! 3 ) * .TRUE • 

GO TO 100 


GRID SIZE CARO 
LSZ»1 

CALL IVAL IICAROilCOLflVEC.LSZ.1002) 

NSMVECtll 

GO TO 100 

ENO CARO 

CALL 1 IVAL 1 1 CARD. I COL ,IVEC,LSZ.ip02) 
FLAGI 51 * . TRUE . 

GO TO 201 


HR TElO.lOU) 
HRITEU6, 10111 
GO TO 999 


CONTINUE 

CHECK IF ALL CONTROL CAROS HAVE BEEN READ 
OD 250 1*1,5 

1FI.N0T. FLAG (1)1 GO TO 321 

GO TO 253 

HR 1 TE 1 16 , 1 01 3 ) 

WRITE! 6,10131 
GO TO 99 9 
CONTINUE 

WRITEll6,ltl 
WRITE! 6, III 
WRITE I 16.1014) 

WRITE! 6,10141 

WR1TEI16.1015) 

WRITE! 6,1015) 

GO TO 680 
WRITE ( 16 , 10101 
WRITE! 6,1010) 


n I I C I D , 1 

50 TO 999 
CONTINUE 


ACAP 


FORTRAN Bl 
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FILE. . . 


AC AP 


FORTRAN Q I 


SUBROUTINE ACUTST 
PURPOSE 

PERFORM THE INITIAL TRANSFORMATIONS OF THE COORDINATES PRIOR 
TO CLASSIFICATION ERROR ESTIMATION 

DESCRIPTION OF PARAMETERS 

A - GOVARl ANE MATRIX ARRAY FOR THE SPECIFIED SUB CHANNE 

AA - DIRECT COPY OF A USED IN INITIALIZATION 

DELTA - ARRAY (IF CELL WIDTHS ALONG FACH FEATURE AXIS 
DET - ARRAY OF DETERMINANTS OF EACH TRANSFORMED COVARIANC 

EV - ARRAV X 0F EIGENVALUES OF THE CURRENT CLASS 

EA - COVARIANCE MATRIX ARRAY FOR ALL THE CHANNELS 

EMA - ML AN VECTOR FUR ALL THE CHANNELS 

HEAD - ARRAY DISPLAYING 1 ITLES 

ir - ARRAY E OETfRMINING THE CYCLIC STRUCTURE OF INOX 
LL - array of THE LCWER COORDINATES of a CELL PU . UUC1 

HA - URIC INAL MEAN VECTORS FOR THE SPEC IF IE C SUB CHANNEL 

MAA - COPY OF MA FOR INITIALIZATION 

NADR - UNHANTbD HF C!)VARUNCE T MATRIX ENTRIES 

• NCC - OES?RrD E CHANNELS l f SUBSET OF T d£ M l9I *}=inRm N A TES 

Ec : 88“.? 5S'ttiHIIMS» c ak u a8II5 

PR - ONE DIMENSIONAL PROBABILITY ARRAY 

0 - WORK VECTOR 

SO - STANOARD F OE VUTlON U ARRAY FOR THE CURRENT CLASS 

UL - UPPER C CORD I NATE OF A SAMPLING CELL 

V - ARRAY COMPRISING PART OF W 

M - THE RESULTING DISCRIMINANT FUNCTION 

HV1 - WCRK VECTOR 

WVZ - WORK VECTOR 

Z PR1 ” I array V ccntaining the coordinates of each sampling 

CELL 

REMARKS 


SUBROUTINE AND FUNCTION SUBPROGRAMS REQUIRED 
GMPRD . ACRE S* PMC »E IGEN.Dt AG.GMTRA 


METHOD 


EACH CLASS IS PROCESSED SEQUENTIALLY. THE ORTHOGONAL 
TrtANSFDRMAl ON THAT DECOUPLES THE FEATURE SPACE IS 
DERIVED BY OBTAINING THE EIGENVECTOR OF THE CL ASS UNDER 
CONC IDER A T 1 ON. THE ENTIRE SET OF S TftT I ST I C S*HE AN VECTORS 
AND COVARIANCE MATRICES* ARE TRANSFORMED RESULTING IN A 
DIAGONAL COVARIANCE MATRIX FOR THE CURRENT CLASS AND GENER 
FORM FOR OTHERS. THE MEAN OF THE NFW COORDINATE SYSTEH IS 
AT THE CURRENT MEAN. THIS PROCESS IS REPEATEO UNTIL THE 
ENTIRE SET OF PUPULATICNS IS EXHAUSTED. 


SUBROUTINE ACUTST U, E y A , A A.S tGMAjEH A , ,MA, MA A, $0 fc OEL T Aj P . PC & V ; WY 
j NADR. NC > AuC * HE AD I 

.♦♦ a **************************** ♦ a ********************************’ 5 " 4 
LOCAL VARIABLES DEFINITION 

****«******* e **** o ****»*********************** < " > ** ,, * < '********** 0 **®® 


FORTRAN BL 


! ritVViSi? 

I NT EGER* A IQINOIMI , I N DX ( NDI MI .N ADR I NS A I . HE AD ( 1 ) . SDIM 

*NCC(SOiMi Ft . , x . • , . ■ / 


INTEGER** FMTI ( 61 /• I 
INTEGER** FHT2(*I / * 2X . ' 
INTEGER** FH T 3 I * I /M2X,',! 


INTEGER** N1 /' 


0'/,N2 /• 


3' ,M 1PE' , U2.A'.' *3X1 •*•»■/ 
0«,'El*.*,'7l'/ 

0* ,'E1A.*, *71 • / 

O'/.Nl /» 0*/ 


C COMMON /CAPCOM/ NO .NO I M,N S , NCLS , NP . I S I ZE , KNT R 

COMMON /ACUCOM/ NTS « N TM * N SS » NS A, NSM » S DI H » NS l ZE * I 5KP 

a® © *********** **** ******* **** ************ *9******** * ***************** 

c COMMON BLOCK VARIABLES DESCRIPTION 

C i SIZE - STORAGE SPACE REQUIRED FOR EACH COVARIANCE MATRIX 

C KNTR - COUNTER FOR THE CURRENT CLASS 

C ND LS - ENnRE C ENTR E tS OF A COVARIANCE MATRIX 

\ r E « ra.si'.aiu ««. 

3!i« : 8 1 8 S!IS ?8'.fiM85 5H 81 ffl »‘i!& 

: 1 8 i!I8 !° SIS8 5?? KISS C & s»c„™ 

l NTS : CORF SPACE ALLOCATED TO COV HAT USING THE SUBCHAN 

C SDIM - TOTAL NO OF CHANNELS AVAILABLE 

C**0*****e*****************»***************' < ' < ’* <l * ,, **' ;,, ‘' ,,,< '' ! ’ < ' ,<l ®*' ,: '**' < ‘ < ‘ <, * ,i ‘ O 
^****,******************************************* < ’ <■*’'' 

70 FORMA T t l H I ( 

71 FORMA T I / I 

72 FORHATI //) 

7* FORMAT! ////I 

536 F0RKAT13SX, 'ANALYTIC CLASSIFICATION ACCURACY PREDICTION - A C A P 

III ^PBRHATilS^SAngSfN&lSRIO CHAgACT|RI|^tg| 0 «^ 

cua FnRMftTl l x. , f»RID CELLS PcK yiNcNi iyij » 

1*1 format! I xJmRANSFORMED FEATURE N SPACE G CHAjActERISTICSO-) 

Kt , 

I** FORMAT! IX. *E IGENVECTORSO' 1 

III FORMAT IX TRANSFORMED MEAN VECTORSO'I c T CLASS 

5*7 FDRMATIZSX.'PRO BAB I L I , | Y OF CORRECT ULA6* 

5*6 ^0RMATI25L'*2***TOtIl|pr6b OF CORRECT CLASSIFICATION****** 

550 FORMAT] IxMlL CLASS STATISTICS HAVE BEEN READ'! 

IG5«NS«*N0IM 

c NI-NUNDIH 

N2 *N2 *SD I M 

Rj 3 Q At 1 X ^ 

IMNSIZE.LT. 51 N3-NS1ZE 

c 

FMTII2l*Nl 
FMT2 1 2 I » N2 
FMT3I2I-N3 


= *#**«»«************************«**** < ‘* <, «*®** 6 * <> ®®***** < ‘*** < ' 0 *‘'**’ > **** 


READ IN THE KEAN VECTORS AND COVARIANCE MATRICES 


£0«»**0*<M 


********************* ********* ******** ♦«*•**•******«**- 


o 

LO 


ORIGINAL PAGE ig 
OF POOR QUALITY. 
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FORTRAN B1 



FILE 


ACAP 


FORTRAN 81 


TPCC»0. 

00 868 KNTR »l.NCLS 


MR I TE ( 6 > 70 1 
WRITE! 6, 721 
C 

C***************** ********* *******9««**********4********6«»*o*9**ea****4 

C 

C RESET THE MEAN VECTORS 

C 

^***0«4>» *0***9 ******* *** «** ******************* 9 ******* *****************4 

c 

00 321 JCLSoltNCLS 
DO 321 KD J K* l * NO I M 

MOSM JCLS-11*N0I H 

MAI HQS+KO I Ml * HA A (HDS+KOlHl 

321 CONTINUE 
C 

c«**** ******** ********** ******* ************** «*««« ******** *************4 

C 

C PERFORM A DECOUPLING TRANSFORMATION ON THE CURRENT CLASS 

C 

C 

C*************** **** ******************************************* ******** 4 

c 


******** »******«««****«« ****«**e*e«04i»*69*«**«*****«««994«*a«** »* 

TRANSLATE THE ORIGIN TO THE MEAN OF THE CURRENT CLASS 
*************************** ********** *0 ******** ******** ****** o * o « 

DO 35 JCLS»1 >NCLS 

IF I JCLS.EQ.KN TR I GO TO 35 

DO 33 KOEN-l.NOlH 

HDS“I JCLS-ll»NDIP 
NOS"IKNTR-I *NDl M 

MA|HOS*KOI Pi=HA( HDS* KD I M I -HA l NDS »KD l N » 


CONTINUE 

CONTINUE 

00 36 KDIH»I,NDIM 

HAI NDS+KO l Ml *0 


4****«e*a*»o«*4*****»ti***** ** ********************** *** ****** *** *4 
FIND THE REQUIRED TRANSFORMATION MATRIX 

»*««**« ************* *4^* ******* ******************** ************* 

FORMAT! I X » A ( 1 PE 1 1 • A* X XI I 

N0S»(KNTR-n*IS!2E*l 

CALL EIGEN { A ! NOS I , R , NDI H ,MV > 

L»0 

DO AO KO IM“1 1 NDI M 


l*l*kdim 

evTkoih) 


KO IMl =A !NDS*L-1 > 


CONTINUE 



FILE. . . ACAP FORTRAN B1 

r tv** ******************************************* ****************** 

C RESET THE CURRENT CLASS COVARIANCE MATRIX 

C ^ tv**** ******************************************** ******** ****** 

c 

DO 901 KDIM-1. (SIZE 

C A<NOS»KDIH-l)«AAlNDS*KDIM-l) 

901 CONTINUE 

C 

C 

C** ******** *********************************** ************************** 

C PERFORH THE TRANSFORMATION ON COVARIANCE MATRICES 

^t************************* *************************** ****************** 

C 

NR l -NOIH 
NCl-NDIM 
NC2-N0IM 
MS A*0 
MSB-l 
ICO*0 
C 

00 701 JCLS-l.NCLS 

C KDS-j JCLS-II*ISI ZE*l 

CALL 01 AG (RjAI*KCSJ»SlGHA(KSS)»XPRIH*HSAtHSBrNRltNCl »NC2 * ICO » 
1 NDIM) 

C 

701 CONTINUE 
C 

£•**«*•** ************************************* **************** ********** 

£ PERFORH THE TRANSFORHATI ON ON MEAN VECTORS 

r 

C****** *************************************** ♦A"************************ 
C 

CALL GHIRA (R,0,NDIM,NDtMI 

NR1-NDIH 

NCl-NUlM 

NC2-1 

DU 702 JCCS-l.NCLS 

CALL *GM PRO M Q^HiwuS I ,M(MDS) ,NRI»NC1 ,NC2I 
C 

£02 CONTINUE 

q* •*••******«******•**••«*• *************************** ****************** 

C PREPEARE THE PRINTED OUTPUT 

£•**«*** ************ *********** *** ************************************** 

C 

KS-2*KNTR 
HR I TE I 6, 5 36 1 

HR I TE I 6, 53 7 ) l HEAD (KS*I> , t-1 ,2 ) 

C HR 1 TE 1 6 1 74 I 

HR1 IE 16,5381 

C HR I TE 1 6* 72 1 

HRITEI6«539) NS 

C HR I TE 16, 72 J 

HRITE(6,540I IGS 

HRI TE 1 6, 74 1 
HR 1 TE 1 6, $41) 


C 
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HR I TE (6, 72) 

HR I TE 1 6, 542) 

HR I TE ( 6 1 721 

HR I TE (6r FHTl ) EV 

HR I TE I 6.74) 

HR ! TE | 6, 544 ) 

HRITE (6,721 
HRITE(6,FMT1) R 

HRITEI6.74I 
HR I TE < 6, 546) 

HR1TE(6,72I 
hriteU.fktii M 

HR I TE (6 , 74 ) 

C* ******************************************************** ****•••*•*••* 

C 

£ START ESTIMATION OF THE PROBABILITY OF M l SCLASS IFIC AT10N 

C********** ****************************************** ****************** 

C 

99 CONTINUE 

CALL PMC IM, SIGMA, P, PC. EV,Z, SO, DELTA, HVt , HV2 ,V, X PRl M ,H, R , 

1 Q.DET ,|NDX, lQ,LL,UL,PR,HEAO) 

C 

C***« *************** ******* ********•******•*•****«*•**»*»*•«•*•••**•**# 

c 

C CONVERT THE RESULTING PCC TO PERCENTAGE 

C******** ****** ******************************************* **•*«*#*•**•• 

C 

PC(KNTH)*100.*PC(KNTR> 

HR I TE 1 16, 72) 

HRITEU6.491 I HEAD I KS *1) , 1 *1 ,2 ) 

49 F0RHAK25X, 'CLASS »,2A4I 
HR i TE 1 16 , 7 1 ) 

HR)TEI16>43) PCIKNTR) 

43 FORMA 7 1 1 OX, • PROBABILITY OF CORRECT CLASSIFICATION- *,F7.3,' - 

l'l 

HRITE 1 16,72) 

HRITEI6,547) PCIKNTR) 

B68 CONTINUE 
C 

£* •*»«•»•**************************•*•****•**«***»•***»*«**«*** ******** 

C FIND THE TOTAL PROBABILITY OF CORRECT CLASSIFICATION 


DO 568 I-l.NCLS 

TPCC-TPCC *PC< 1 1 
CONTINUE 

TPCC- T PCC /FLOAT (NCLS) 
HR I TE ( 6 , 74) 

HR I TE i 6, 546) TPCC 
HRITE 16, 721 


CONTINUE 


RETURN 

END 


SUBROUTINE ACRES 


l — ' 
O 
■C- 
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LOWER COORD. OF A CELL 

T R ANS FOR MED MEAN VECTORS 

GRID POINT VECTOR 

CL AS S I F I CAT I OH ACCURACY RESULTS 

PRCBAB I L I T Y ASSOCIATED WITH EACH GRID CELL 

WURK VECTOR 

E I CENVCC TORS 

CLASS STANDARD OEVIATION VECTOR 
PART OF IHE DISCRIMINANT FUNCTIONS 
WORK VECTOR 
WORK VECTOR 
0 1 SCR I M [ NAT FUNCTION 


XPRIM- WORK VECTOR 


- COCRD OF CELL CENTERS 


REMARKS 


SUBROUTINE ANO FUNCTION SUBPROGRAMS REQUIRED A 

A 

01 AG. GRPRQB, MI NV A 

METHOD A 

A 

THE PRIMARY INPUT TO THIS SUBROUTINE IS THE TRANSFORMED A 

MEAN AND COVARIANCE MATRICES. THE STRUCTURE OF THE SAMPLING A 
GRID IS OETERMIIUD BY COMPUTING THE CELL WIDTH ALONG EACH A 
DIMENSION ANO THfc COORD OF THE CENTER OF EVERY CELL THKDUGHQUA 
THE GRID. THE ENURE GRID IS SCANNED AND NCLS OISCRIMINAT A 
FUNCTION IS CALCULATED FOR EACH CELL. USING MAXIMUM LIKELIH A 
OJO RULE A CELL IS ASSIGNED TO EITHER THE CURRENT CLASS OR A 
OUT 5 1 DE ( NE ED NOT KNOW EXACTLY WHICH CLASS). SUBROUTINE A 

GRPROB CALCULATES THE HYPERVDIUMC JNDER THE PDF AND OVER THE A 
GRID CELL. AFTER ALL THE POINTS ARE EXHAUSTIVELY TESTED THE A 
ELEMENTARY UNITS OF PROBABILITY ARE SUMMED AND STORED IN THE A 
PC ARRAY AND RETURNED TO THE CALLING PROGRAM. A 


C A 

C A 

SUBROUTINE PMC <M, S IGHA , P , PC ,E V , Z , SO, DEL T A , WV 1 , W V2 , V , X PR 1 M,W , R , A 

1 Q.DET ,INDX, 1 C.LL.UL, PR, HEAD) A 

C A 

C A 

a*#* e*u* o .**«* e •» e*** * o, ?■* «»*«*«**«******* e****** * a 

C A 

C LOCAL VARIABLES DEFINITION A 

C A 

C A 

REAL MIH.LLU I ,ULUI ,PR| 11 , A 


I'm n< WtLI. I W lULI II ll’l'lll t 

1 SIGHAIll.ZlNDIM.NSI.SUt II. DELTA 111, 

2 Pll I ,WC 1 1 1 ,WV1 1 h.WVZC 1 I , PClTS , VI li • 

3 XPR 1 HINDI Ml, WI1 I ,R(l l,Qli» ,DET( 1 ),EVI II 

INTEGER* A GS, HE AO 111, IQ 1 1 I , I NDX 1 1 ) 

COMMON /CAPCOH/ ND ,ND IM,N S, NCLS , NP, I S I ZE , KNT R 


GS“NS**NDIM 

ICD«1 


C* ««« «*«*«* at t> «****• *»**«*» ?#****»* *•*****•*»***«*•*« A 

C A 

£ FIND THE VARIANCES ALONG EACH FEATURE AXIS A 

C*»**o**o »**<■** **<>»*» **«**» ♦*****«a*»CM,*»*»* 90** «*** **«*»*«* *«* » + a 

C A 

DO 706 KDIM*1,N0|M A’ 

C A 

SO ( KO INI “SORT IEVIKOIMH A 


CONTINUE 
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C******* ****** ************* *********************************** ********** 

C INVERT THE TRANSFORMED COVARIANCE MATRICES 

C**«***A* ♦*•«*• ******************************************* ************** 

° DO 70 T JCLS-l, NCLS 

CALL ' HINV ! SI GMA I KDS ) ,NCIM»OET I JCLS I » WV1 ,WV 21 
707 CONTINUE 

33 FDRHATI1X,3I 1PE12.A.2XJ1 

32 FORMAT! IX, 2( 1PE12.AI I 

C 

C**«* A************************************************ ****************** 

C SAMPLE THE FEATURE SPACE BY A BINOMIAL APPROX. TO NORMAL O.F. 

C ♦**•***♦ ************************************************************** 

C 

S 

c *************************************** a******* ****************** 

C FINO THE WIOTH OF A SAMPLING CELL ALONG EACH DIMENSION 

£ ***************************************************************** 

c 00 AA KD I M*1 * NDI M 

C DELTAIKDIM)«2.*SD(KDIM) /SORT I FLOAT I NPI I 

AA CONTINUE 

C 

Q *************************************************** ************** 

C FIND THE COORDINATES OF EACH AND EVERY SAMPLING CELL CENTER 

£ **************** ******** *********************** **** ************** 


00 30 KD I M-l »ND IM 
DO 30 I-l.NS 

C ZIKOIM,II-2.*SO(KDIMI*!FLOAT!l-n-NP/2)/SQRT(FLOAT(NPn 

C 

30 CONTINUE 

C 

£*»*•«•*• A***************** ******************* *♦******♦♦♦**♦*♦* ********* 

C* CALCULATE THE CUADRATIC OISCRIHINANT FUNCTIONS * 

^^«***** ******* ******** ******* *********************** **** ************** 

° ICNT-0 

00 AO A KDIM-l.NDIM 

AOA IN0X(K0IM)»1 

C 

5 ***************************************************************** 

C PERFORM A ONE-TO-N DIMENSIONAL MAPPING OF POSITIVE INTEGERS 

r a**************************************************************** 

DO 520 I * 1 »N0 1 M 

c IQI n»NS**(N0IM-n 

c 

520 CONTINUE 

DO UO J-I.GS 

IFIJ.EQ.1I GO TO 30 


C 
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DO 122 KDIM-1,NDIH 

C IFIKODIJ-I.IQIKDIMH .EQ.Ot INDX (KOI M I - 1 NOX ( KDI N I *1 

IF! INOX(KDiH) .GT.NS) INDX(K0IH>*1 
C 

122 CONTINUE 
3« CONTINUE 
C 
C 

c *«••*****«*****•*•*«•«**•*••*****•*******•••*••«*•«*••**•••*****< 

H COMPLETE THE DISCRIMINANT FUNCTION CALCULATION 

£ •*•***• ******************** ********** ****** *•***«*•*•**«**••*•*•< 

00 12A JCLS-I.NCLS 
DO 126 KDIH»1,NCIM 

MDS*( JCLS-l )*NDl M 

P t MOS*KO I HI-2 IKO IM, I NDX ( KD I H»> -Ml HDS *KD IM1 

126 CONTINUE 
12 A CONTINUE 

C 

00 12B JCLS-l, NCLS 
C 

HDS-! JCLS-ll*NOlM*l 
KOS-I JCLS-l 1 *ND+ l 
MSA-0 

CALL°DIAG (PIHOS t .SIGMA IKDSI ,V ( JCLS I ,XPRIM, MSA, HSB.NR1, 

I NC I , NC2 , 1 CD, ND I Ml 

' HIJCLSI-VC JCLSI + ALOG!OETI JCLS1 I 

C 

128 CONTINUE 

^*** ************************************************* ***** ************** 

C* CHECK THE CONDI 1 1 ON FOR CORRECT CLASSIFICATION 

£**«« ************************************************* ******** •»*•***•** 

C 

TE1M-1.0E 10 
00 I3A JCLS-l, NCLS 

IF<W( JCLS) .LE. TERM! TERM-WIJCLS) 

C 

1 3 A CONTINUE 


£•*****# $**•***♦ *«****?**** a<r******##*o****$**#**#***<ii><«<i*»***e »*•**♦*** 

£ FIND THE ELEMENTARY UNIT OF PROBABILITY 

r 

L* ************************ ************************ ******** ************** 

L IF! TERH.NE.WIKNTH)) GO TO 120 

ICNT«ICNMl Rpft o B , 1 nox , OE lta,SO,NS,NDIM,Z,PR,LL,UL) 


CONTINUE 

PRP-l. 

DO 710 KD IH-l « NCI H 

PRP-PRP*PR!INDX(K0IM)1 

CONTINUE 

PCI KN TR 1 • PC IKNTR I *PRP 

CONTINUE 

RETURN 

ENO 


C 


o 
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SOB ROUT 1 N GRP RUB 
PURPOSE 

COMPUTE THE GRID PROBABILITY 

DESCRIPTION OF PARAMETERS 

DELTA - GRID CELL DIMENSIONS 

INDX - POINTER ARRAY 

LL - LOWER CHORD OF A CELL 

NDIH - NO CF SELECTED SUBSET CHANNELS 

NS - NO cr CELL PER DIMENSION 

PR - PROBABILITY ARRAY AROUND EACH CELL 

SD - STANDARD CEVIATICN VECTOR 

UL - UPPER COORD OF A CELL 

l - CELL CENTER COORD ARRAY 


REMARKS 


SUBROUTINE AND FUNCTION SUBPROGRAMS REQUIRED 


HETH00 

THE NORMALIZED CLASS DENSITV FUNCTION IS INTEGRATED OVER 
A HYPERVOLUME WHOSE SIZE IS DETERMINED BY 'DELTA*. 


SUBROUTINE GRPRCB 1 1 NDX.OELTA, SO, NS ,NDIM, Z ,PR,LL ,UU 
REALPWINOIM.NSj.OEL TM1 I ,PRU> .SOI U.LU 1 > ,UL< I) 

00 100 KD!M»l,NDiH 

ULIKDIMI=Z(KDIH, INOX IKOIMJ I ♦DELTA (KOI HI ZZ. 

LLIKOIHI-ZIKDIM, INDXIKOIMI 1 -DELT AT KD I MI /Z . 

UL(K0IH»»UL|KOIH I/SD IKOIHI 
llIkdimt«llIkoim izsdIkdihi 

i PR j j NDX (KOI Ml ) *0 »5«( ERF (UL I KDI MI/SQRT 12 .1 1- ERF I LL( KOIM IZSQRT 1 

CONTINUE 

RETURN 

END 


SUBROUTINE OIAG 
PURPOSE 

TO PERFORM THE FOLLOWING OPERATION. I QTRANSPO SE I * X* 1 0 1 
USAGE 

CALL OIAG IQ, X, SIGMA, XPRIM.NSA, MSB, NR1, NCI, NC2,ICDtN> 

DESCRIPTION Or PARAPFTERS 
0 -THE FIRST MATKI X 

X -THE SECOND MATRIX 

SIGMA -THF OUTPUT MATRIX 

MSA -FLAG FOR THE STORAGE MODE OF FIRST MATRIX 

MSB -FLAG FOR THE STORAGE MODE OF SECOND MATRIX 

NRl -NO DF ROWS IN THE FIRST MATRIX 
NCI -NO OF COLUMNS IN THE FIRST MATRIX 

NCZ -NO OF COLUMNS IN THE SECOND MATRIX 

ICO -FLAG TO RESET NR AND NC PARAMETERS 

H -NO Of DIMENSIONS 
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RETURN 
ENO 


SUBROUTINE TPRO 

PURPOSE 

TRANSPOSE A MATRIX AND POSTHULTIPLY BY ANOTHER TO FORM 
A RESULTANT HATRIX 

•USAGE 

CALL TPRD(A,B,R,N,H,MSA,HSB,L) 

DESCRIPTION OF PARAMETERS 

A - NAME CF FIRST INPUT MATRIX 

B - NAME OF SECOND INPUT MATRIX 

R - NAME OF OUTPUT MATRIX 

N - NUMBER OF RCWS IN A AND 0 

M - NUMBER OF COLUMNS IN A ANO ROWS IN R 

HSA - UNfc DIGIT NUMBER FUR STORAGE MODE OF MATRIX A 

0 - GENERAL 

1 - SYMMETRIC 

2 - DIAGONAL 

HSB - SAME AS MSA EXCEPT FOR MATRIX B 
L - NUMBER OF COLUMNS IN B ANO R 

DC U AD K S 

MATRIX R CANNOT BE IN THE SAME LOCATION AS MATRICES A OR B 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
LOC 

METHOD 

MATRIX TRANSPOSE OF A IS NOT ACTUALLY CALCULATED. INSTEAD, 
ELEMENTS IN MATRIX A ARE TAKEN COLUMNWISE RATHER THAN 
ROWWISE FCK MULTIPLICATION BY MATRIX B. „ „ _ 

THE FOLLOWING TABLE SHOWS THE STORAGE HOOE OF THE OUTPUT 
MATRIX FOR ALL COMBINATIONS OF INPUT MATRICES 

A B ^ 

GENERAL GENERAL GENERAL 

GENERAL SYMMETRIC GENERAL 

GENERAL DIAGONAL GENERAL 

SYMMETRIC GENERAL GENERAL 

SYMMI.TRIC SYMMETRIC GENERAL 

SYMMETRIC DIAGONAL GENERAL 

DIAGONAL GENFR4L GENFRAL 

DIAGCNAL SYMMETRIC GENERAL 

DIAGONAL DIAGONAL DIAGONAL 


CALCULATED. INSTEAD, 
ISE RATHER THAN 

MODE OF THE OUTPUT 
ATRICES 

R 

GENERAL 

GENFRAL 

GENFRAL 

GENERAL 

GENERAL 

GENERAL 

GENFRAL 

GENERAL 

DIAGONAL 


SUBROUTINE TPRDU.B.R ,N, M , MS A, HSB , L 1 
DIMENSION All) ,e(l),R(li 

SPECIAL CASE FOR DIAGONAL BY DIAGONAL 

HS*HSA*1 0*NSB 
IFIKS-22I 30,10,30 
10 DO 20 1 * j ,N 
20 R I I )>A( I I *U I 1 1 
RE TURN 

MULTI PLY- TRANSPOSE DF A BY B - 

30 IR«1 

DO 90 K* 1 «L 
DO 90 J*1.M 
Rt | R ) >0, 0 
DO BO I » 1 . N 
IF I MS I 40,60, AO 


iMnu *»u,ou»^u 
40 CALL L0CII,J,IA,N,M,HSA 
CALL LOCI I , K , I B, N, L »H SB , 
IF I T A I EO, 60,50 
50 IF i 18 I 70.60,70 


5U (MIDI , Ul HU, /U 

60 I A >N» I J“ 1 1 ♦! 

1 8 *N* IK-t )♦! 

70 RtlRJ-RI IR)*AUA|O01!0I 
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80 CONTINUE 
90 IR-IR*1 
RETURN 
ENO 
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£**04** *04$ 4 Q* OOOOOO 0**0 *00000$ *** $**$$*$*0**$ *0*00*0**000 0***0 ♦OOfr***'* 

C SUPERVISOR FOR THE STRATI F I |0^ POSTER ! OR CLASSIFICATION ACCURACY 
C WRITTEN 01/29/7B 0. J. W1ERSHA 

C • o 40 *»»***♦ **»«®»SA******I>4>*« ***»® ****** *** ******** ***********+* 

£ AMENDEO 10/31/78 B. G. KOBASSERI 

C DESCRIPTION AND PURPOSE 

C ' THE STRATIFIED POSTERIOR ESTIMATOR - SPEST- IS SIMILAR TO TH 

C ACAP PROCESSOR IN PURPOSE. IT PROVIDES A SET, OF CUSSIFICATIO 

r ACCURACY ESTIMATES FOR A BAYES CLASSIFIER WHEN THE POPULATION 

C STATISTICS ARE MULTIVARIATE NORMAL. ALTHOUGH THE ALGORITHM 

C EMPLOYS THE LI KbLIHOOO PRINCIPAL. IT DIFFERS FROM ACAP APPR0C 

C BY USING A RANDOMLY GENERATED OATA BASE. FURTHER DETAILS 0N„, 

(j THE THEORY CAN BE .FOUNO IN THE LARS FINAL REPORT.NOV 30 ,1978 

C DESCRIPTION OF CONTROL CARDS 

C PSPEST 

C THIS CARDS SPECIFIES THE PARTICULAR PROCESSOR REQUESTED 


CHANNELS 


THE DESIRED SUBSET OF THE AVAILABLE CHANNELS IS GIVEN HE 
I T IS - IMPORTANT TO REMEMBER THAT THE NUMBERS APPEARING 
ON THIS CARO IS THE ORDER OF THE SELECTED CHANNELS NOT 
THEIR ACTUAL NUMBER. FOREX AMPLE ,, I F THE AVAILABLE 
CHANNELS ARE 8,9,12,19 AND CHANNELS 8,9 AND 19 ARE REQUE 
THIS ‘ C ARO SHOULD READ 1,2,9. 

CLASSES 

THIS CARD SPECIFIES THE NAME OF EACH CLASS. EACH NAME 
MUST BE PLACEO IN A FIELD 7 CHARACTERS LONG FOLLOWED 
BY A BLANK. THE CONTINUATION CARO , IF REQUESTED, MUST. ) 
VE THE NAME 'CLASSES* IN THE BEGINNIG FOLLOWED BY THE 
REST Of THE NAMES. 

END 

SIGNALS THE END OF THE CONTROL CARDS 
*«• STAT D.ECK FOLLOWS IMMEDIATELY. MUST BE IN CHARACTER FORMAT 
REHARKS 

THIS PROGRAM IS CURRENTLY CAPABLE OF PROCESSING UP TO 20 
CLASSES AND B. SPECTRAL BANDS; 

HOW TO RUN THE, PROGRAM 

TYPE 'GETDISK OHSYS* TO ESTABLISH THE PROPER LINKS. THE RET 
ER-'HLE CONIAI NS ONE DECK CONSISTING OF THE CONTROL CARO 
AND LARSYS STATISTICS DECK. TYPE * SPEST' FOR THE PROGRAH; 
EXECUTION. 

EXAMPLE OF A CONTROL CARD SET UP 


* SPEST 

CHANNELS 1,2,9 
CLASSES BARE SO I CORN 
END 


PASTURE WHEAT 


C® ******* »♦*»♦ «•***♦ *«*«»***•***♦ ***«♦»*♦*<■**$ * ****** *****************> 


RE AL*9 
REAL*9 
R£AL»9 
RE AL*9 


S SS ! IS i :E8l } S t?8i stl J ( 751 f I IS! ;|S! i2f i Kio . 
siiskHj}u c x88S!n:Sb5 ,ll5 “ , H Sil,s< " 
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5 ?l i? l / ^ is PE?i *CC AS • . • END* V 

1 ' :ra}!' 1 C ^^RD/' LARS'/ 

INTEGER * 2 ICSEL ( 30 ) ,NCj JO ) ,NCC I 3 0 ) 

COMHON L /SPCaH/ , NTHtNTstililE,NSI 2 E,NSA,SDIM,NDlH 

C 

1010 FORMA T I 1 X ! '£ RROR IN CTLWRO. EXECUTION TERMINATED.') 

10 2 FORMAT IX 'ERROR IN CMANElI EXECUTION TERMINATED') 

1013 FORMAT! IX, 'HISSING CONTROL WORD. EXECUTION TERMINATED') 

1019 FORMAKIX, 'SPEST PROCESSING STARTED') 

**• 

c 

c oecooe control cards 

C******* *********************** 6*********************************** 0 ** 

c 

DO 777 1 * 1,30 
ICSEL 1 1 ) ”0 
NCI I 1*0 

777 CONTINUE 
C 

DO 778 I * 1,90 

778 ICSETI 1 1 *- 50000.0 
C 

100 CONTINUE 
LS 2*9 
l£K *0 
INRO* 5 

C CALL CTLWRO II C ARO , I COL.L l ST ,1 S2, 1C0DE , I NRD, IER ) 

IFUER.NE.O) GO TO 1001 
GO TO 199 , 101 , 102,103 ), ICODE 
99 FLAGU )“. TRUE. 

GO TO 100 
C 

C CHANNELS CARO 

101 CALL CHANEL UCARD, ICOL, NCR , ICSEL, ICSET ,NCC ,900 ) 

FLAG I 2 I* .TRUE . 

NO I M« NCR 
GO TO 100 

E CLASS NAMES CARO 

C 

102 00 10 1*1,20 

10 HEAD! II-ICARDUI 

FLAG! 31 *. TRUE. 

GO TO 100 

E ENO CARO 

103 FLAG! 9 ) -.TRUE.' ' 

GO TO 201 ' ■' 


CONTINUE 

CHECK IF ALL CONTROL CARO 5 HAVE BEEN READ 

m ? not! flJgi i i ) go to 32 i 

WR!IEU6?l013j 
WRITE! 6,10131 
GO TO 99 9 
CONTINUE 

WR I TE 1 16, 11 ) 

WRITE! 6,11) 

WRITE 16, 1C19I 
WRITEl 6,1019) 

WAITE 1 16 , 1015 ) 


O 
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CALL MCOVPlNtAK(l«ll« CO VI 1 « I )) 
Eall'spIstmim.n.phi.p.am.cov.pr.pc.op.cov^ 

1 Y.TEI.PEL.CDVU.PX.SDET.NADR.NC * NCC . E X AM . EXCO V . 

2 HEADI 


PC« 100 .*PC 

60 F 0 RHAT(J// 2 / 30 X, 'OVERALL PROBABILITY OF CORRECT RECOGNITION - *,F 1 

1.31 

999 STOP 
END 

C - * 

C SPEST IS AN ESTIMATOR OF THE CLASSIFICATION PERFORMANCE FROM A 
£ GIVEN SET OF STATISTICS FROM M CLASSES. , THE ESTIMATOR IS A 
C STRATIFIED POSTERIOR ESTIMATOR (REF. WH 1 TSITT ANO LANDGREBE). 

C THE PROBABILITY DISTRIBUTIONS APE ASSUMED TO BF MULTIVARIATE 
C GAUSSIAN 
C 

C 19 JANUARY, 197 B 


^SUBROUTINE SPE S™ « .N^PH I jP , AM .COV , PRjPC J « g 6 g 0 VT .CAM . DE T ft C0 V I N . . 

2 HEAT) 

INTEGER** NADRINSAl.SClM 

INTEGER** FHTllAI /M2X,', 2 * < 

INTEGER** FMT 2 I* I / , ( 2 X,‘» I 0 1 , ' E 1 * . • . 1 7 I ' / 


rpal***opimKpimj%r MUEXAH?S 0 IM,M 1 .CXCOV INSIZE .MI.COVTI IS 1 2 E I 
REAL** GAM (h!h } , PH 1 1 N <N,M I .OCT I M 1 ,COV IN l ISlZE.Hi 
REAL** YIN) .TEliN) .nELINI .COyUIN.N) 

saw KlSiJiiS SM 


LIST OF VARIABLES 

K ■ NUMBER OF CLASSES 

N * NUMBER CF D I PENSIONS , 

Pill « APRIORI PROBABILITIES CF CLASS I 
PRIU « CLASS CONDITIONAL PERFORMANCE 
PC - OVERALL PERFORMANCE , 

COVl’i 1 ,}) ’ CUVARIANCE MATRIX OF JlaSS IISTOREO IN UPPER TRIANGULAR 
FORM J 


CO oooo* 0 *o* ******** *♦**♦*♦♦♦*♦*♦♦*♦**** ••♦***** 0 * 000000**0 oo ooo ooooo • * 

C READ IN THE MEAN VECTORS AND COVARIANCE MATRICES 

C ******* ******♦*«*♦******♦***♦*♦*'***♦******•******•*** ****♦*•*******'*•** 

C Nl*Nl*S01H 

N 2 *N 2 ♦ 5 

IF(NSUE.LT*5I N2-N2*NSUE 

° FMT It 2 1 *N1 

FH T2I 2 1 *N2 
C 

READ°(i*FHTll IEXAHIJ.II, J-l.SDIHI 
30 CONTINUE 

IlEAD 1 I5 ,£Ai 2I lEXCOVIX.IJ .K-l.NSIZE) 


KE»U IJtC 

CONTINUE 


H 1 

O 
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coo***** **»*a* «*♦»*» **0»*e* *o»*****»*si«** ****** ******* ****************** 

C 

jj SELECT THE REOUIRGO SUBSET OF THE CHANNELS 

Coco***** *************** *8 *0*** ******************* ********************** 

c 

DO 803 JCLS«1,M 

c 

K”0 

DO 80Z K01K-1,SDIN 
DO 801 J«1,N 

IfSkBihIIqI NCC I j ! { AMI*! JCLS)»EXAN(VCCI Jl, JCLS) 

C 

801 CONTINUE 

802 CONTINUE 

803 CONUNUE 

C 

C**«* ****** *•*««***« *********** ****** ************* ************* ********* 

C OBTAIN THE ADDRESS OF THE SELECTED ENTRIES INTO COVARIANCE MATRIX 

Ce*e*«***««*a«»*»**9* ********** *************************** ***#■»****♦***' 

C CALL ADRES I NAOR.NC, SDI H, NDI N I 
C 

K-0 

DO 603 JCLS»1,N 


DO 602 I 
DO 601 J 


-ItNSIZE 

•liNSIU 


IFU .EQ.NADRI JU GO TO 602 
CONTINUE 

COV<X,JCLSI*EXCOVt 1. JCLSI 


CONTINUE 

CONTINUE 


IX » 9AT913 
NCT * N* (N»U F? 

COMPUTE EIGENVALUES AND EIGENVECTORS FOR EACH MATRIX 

MV ■ 0 
EPS • l.Ot-6 
00 100 IJ*l,H- 
DO 56 ]«1,NCT 

CDVTIll a COVI I * I J I , , „ „ „„„ 

CALL EIGENICOVT.PHHl ,1,1 Jl.NlHVI* 

L » 0 

00 60 I * 1 ,N 
L * L ♦ I 

GAM 1 1 ,IJI » COVTILI 

COMPUTE DETERMINANT AND INVERSE OF EACH MATRIX 
DO 65 I-l.NCT, 

covim ■ coyii.iJt .... „„„ 


CALL SMI NV ICOVT ,N.DET IIJI ,MV .EPS.IERI 

I F I 1ER) 1000, 70, lOOO 

CONTINUE 

SDETIIJI ■ SORT IDE T ( I J> ) 

DO T5 lal.NCT . 

COVINll.lJl * COVTIII 
CONTINUE 
MV ” 0 

DO 105 I-1.H 


SORT IOE T 1 1 J> 1 


COVTIII 


DO 105 IM.H 
105 QPIII a 0.0 
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LOOP ON CLASS ICL 

PC ■ 0.0 
DO 500 ICL-l.M 
AVEO - 0.0 

LOOP. ON THE NUMBER OF SAMPLES 

NS « 1000 
DO 300 I J»l,N$ 

GENERATE Y VECTOR FROM CLASS ICL 

00 110 I-l.N 

CALL RANOUIIX,lY,XPI 

IX » IT 

CALL NDTRI1XP,Y(II,XD,IERI 
0 CONTINUE 

CQHPUTE CONDITIONAL PROBABILITIES FOR EACH CLASS 

00 200 JCL»l ,M 

1 F I JC L .EQ. ICL I GO TO 180 
00 130 I * 1 ,N 

TE1I I I « 0.0 

DELI] ) » AMI I* I CL I - AMII.JCU 
OD 130 J»1,N 

T E 1 < I ) * TEl 1 1 I a SQRTIGAHI J.ICLI miJI*PHII I,J, ICL I 
0 CONTINUE 
JJ * 0 

00 1 A 0 1*1, N 
DO I AO J»l,I 
JJ « JJ * 1 

COVUI I , J I • COVIN! JJ, JCL 
CDVU(Jtl) - CDVINI JJ. JCL I 
0 CONTINUE 
Z1 « 0.0 
Z2 a 0.0 
13 » 0.0 
DO 150 I • 1 ,N 
DO 150 J»1,N 

Zl * Z1 - 0.5*TE1(I I *COVU II,J)*TE1IJI 


Zl • Z1 - 0.5*TE1II I *COVU (I,J)*TE1(J) 
Z2 a Z2 - Teim»CUVU(I,J)*DELtJI 
Z3 « Z3 - 0.5+DELll) ®COVU (1,J)*DEL(J) 


continue 

Z SUM a Zl ♦ Z2 ♦ Z3 
IFIZSUM .LT. -1001 GO TO 190 
BETA " P t JCL 1*1.0 

PXIJCL) • BETA'*DEXP(Zl*Z2*Z3)/SOETIJCLI 

IRIPXIJCL) .EQ. O.Ot WRIT El 16, 9191 ICL, JCL.ZSUM, SOET ( JCL l ,PXI JCL 

CONTINUE 

GO TO 200 

CONTINUE 

ZO » 0.0 

DO 185 l'«l,N 

ZO « ZO- 0.5*YI I l*Yt l I 

CONTINUE. 

IF I ZO .LT. -1001 GO TO 190 

BETA a p( JCL 1*110 

PXIJCL! • BETA*OEXP(Z0l/SDETIJCLI 

1FIPXIJCLI .EQ. 0.01 WRITE(16,919) 1 CL, JCL ,Z0, SD ET I JCLI , PXC JCL I 
GO TO 200- 
PXIJCL) * 0.0 
CONTINUE 

FORMAT I 5X, 215,3612.41 


CHOOSE THE LARGEST 


BIG » -1000 
DO 220 I»1,H 
IPIPXIti .GT. 


BIGl LOC a 


1FIPXI I t .GT. BIG) BIG a 

Continue 


continue 
oen - L o.o 

OD 230 t*l.R 
DEN DEN ♦ PXI II 


220 
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CONTINUE 
Q * BIC/OEN 


QPUOCI - QPILOC! ♦ PUCLI*Q/PILOCI 

CONTINUE 

CONTINUE 

HRITEI6. Ill 

FORMAT! ^5X, •STRATIFIED POSTERIOR ERROR ESTIMATOR 1 ! 

FORMA T 1 1HI I 

DO 510 t CL * 1 f M 

PR! 1C L 1 * QPl ICU/FLOATINS! 

PC « PC * p< ICL l*PR( ICU 
WR ( TE ( 6 » 72 1 
FORMA T l / / / I 
DiPL*2*ICL 

WR I TE 16, 53 7) (HEAD (OS PL*- 1 ) , I -1,21 
F0RKATIA9X, ‘CLASS '.2AA! 

PC PR* 100 . * PR ( I CL) 

WR ITE lb. B2I 
FORHATI//) 

WRITE! 6, AB2I PCPR _ 

FDRMA T 1 2 5X . ' P ROBABI L ITT OF CORR! 
LS1FICATI0N* ' ,F7.3»' 1 1 ! 

CONTINUE 

RETURN 

WR ! TE < 6. 1100 1 1 ER 

FORMA TtlOX,**** INVERSION ERROR 1 1 ,12, '>***' I 

RE TURN 

END 


CORRECT CLAS 


SUBROUTINE ADRES 
PURPOSE 

FIND THE DESIRED SUBSET OF A COVARIANCE MATRIX 

DESCRIPTION OF PARAMETERS 

NADR - ADDRESS ARRAY OF THE ENTRIES TO BE DELETED 
NC - ARRAY OF THE CHANNELS TO BE DELETED 
SOIM - TOTAL NO OF BANDS SUPLIED 
NDIM - OESIRED SUBSET OF SDIM 

REMARKS 


SUBROUTINE AND FUNCTION SUBPROGRAMS REQUIRED 


METHOD 

USING A ONE DIMENSIONAL STORAGE MODE FOR £*CH COVARIANCE 
MATRIX IN AN UPPER TRIANGULAR FORM, THE PARAMETRIC ADDRESStS 
FOR EACH ENTRY IS OERIVEU ANO USING THE NC ARRAY THE LOCATION 
OF ALL THE ENTRIES THAT LIE IN THE UNWANTED LINES AND COLUMNS 
IS COMPUTED AND STORED IN NADR. 


SUBROUTINE ADRE SINAOR ,NC, SOI M.NDIM) 

INTEGER** NADR 1 1 ! • POI M»SO 1H 
INI EGER*2 NCU! 

PD I M» SD I M-ND I H 

IF tPDIM.LE.O) GO TO ICO 
00 20 J-l.PDIM 
N1«NC(JI 
DO 10 l»l,Nl 


. . SPEST FORTRAN Bi 

NADRIK)* (NCI J)*INCIJ)*1)/2)*I 

CONTINUE 

CONTINUE 


DO *5 J-l.PDIM, 
NF I N» SDl M-NC ( J I 
DO 35 l*X »NF IN 


NACRIK)«MNCIJ!*I-X)*«NCI Jl*l-2I/2»*NC(J) 

CONTINUE 

CONTINUE 

RETURN 

END 


p P 
^ 1 # 

cj3 

gg 

§ g 

S H 

% o5 


a. a. a. a. a. o. a. cl cl q. & tx. cl cl cx a. 
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r, # *0«(,*(M«*««(.»«*»«t4H«»<lM»*iH*«««*««***W****W******” <,M *'* 

C SUPERVISOR FOR A TWO DIMENSIONAL SPATIAL CORRELATION 
c MRIUEN 10/0A/78 Bl JAN G. HOBASSERI 

ce**0'>*o*******««*«**P**< , **' > ***** ,, ' <,,M ** ,> * <,,< ’‘ > *** t " , ‘ l! '***'*********“ > *** ,> * , ' < ' 

C 

r DESCRIPTION AND PURPOSE 

r r H ,c punr.RAH IS A TWO DIMENSIONAL SPATIAL CORRELATOR THE 

r PRIMARY OUTPUT OF VHICH IS A NORMALIZED SPATIAL CORRELATION 

C MATRIX FOR ANY SPECIFIED AREA. 

r THF USER SPECIFIES THE COORDINATES OF HIS DESIRED AREA IN THE 

r FORM OF INITIAL AND FINAL LINES AND COLUMNS ALONG WITH "THE 

SS^irfAii itft T aSS S&sBA84*NASBK e N#K°irattH«o: u, * TlOM 

\ ISfchSg&?5Ai H Pif s SifiW°^^iWgKi8t R SV£J l RiT H ffl R 18^NB|Ti*LLv 

r DROPPING FUNCTION TO THE EXPERIMENTAL DATA USING A WEIGHTED 

£ LINEAR LEAST SQUARES FIT TECHNIQUE. I n^cc cn^ ^ 

C ASSUMPTION IS NOT VALID, THIS OPERATION IS BYPASSED. 

r THF PURPOSE HERE IS TO DEVELOP A HARKOV MODEL FOR THE 

l SPATIAL CORRELATION FUNCTION OF THE MSS OATA. 


DESCRIPTION of CONTROL CARDS 
♦CORRELATE 

THIS CARD SPECIFIES THE PARTICULAR PROCESSOR REQUESTED 
INPUT RUNI.) , TAPE! . 1»F1 LEI .» 

THE INPUT RUNTABLE FROM WHICH DATA IS READ 
BLOCK LINE!. «.) .COLUMN! .,. ) 

SEGMENT TO BE CORRELATED 
FUNCTION 

EITHER AUTC OR CROSS FUNCTION CAN BE SPECIFIED 
CHANNELS 

SPECTRAL BANOS USEO IN CORRELATING 
SAHPLELAG 

CROSS TRACK LAG USED IN ESTIMATING THE .CORRELATION FUNCT ION 

EXPRESSIONS A PERCENTAGE. OF TOTAL NO OF SAMPLES IN 


THE AREA. THI S CARD 


OPTIONAL. 


LINELAG 

SAME AS SAHPLELAG EXCEPT FOR LINES (ALONG TRACK I 
EXP OF IT 

THIS CARO STARTS THE EXPONENTIAL FITTING PROCEDURE. OPTIONAt 


END OF CONTROL CARDS 
REMARKS 

this p?Si5t M i.iSoi u 5S?SIEfe v AEk p TBfe e s85R6B?1J511 ^Se^vhS^callv 

• SImENsIoNEO. ANY ENLARGEMENTS CAN BE ACCOMPLISHED BY ALTER. N 
THE DIMENSIONS OF THE ARRAYS IN THE MAIN PROGRAM. 
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HOW TO RUN THE PROGRAM 

USE THE COMMAND 'GETOISK DHSYS' TO ESTABLISH THE PROPER LINK 
THEN TYPE ^CORRELATE. 

EXAMPLE OF CONTROL CARO SETUP 

♦ CORRELATE 

INPUT RUN! 7A02B500) . TAPE 12689) .F ILEI 3) 

BLOCK LINE! 1,25), COLUMN!!, 25) 

FUNCTION AUTO 
CHANNELS 2 
SAHPLELAG 25 
LINELAG 25 
EXPOFIT 
END 

OR FOR CROSSCORRELATION. 

♦CORRELATE 

INPUT RUN I 7402 8500, TAPE 12689 I.F ILEI 3) 

BLOCK LINE 1 1»ZSI * COLUMNIl, 25) 

FUNCTION C R CSS 
CHANNELS I, 3 
• SAHPLELAG 25 
LINELAG 25 
EXPUFIT 
ENO 

IF SAHPLELAG AND LINELAG ARE LEFT OUT. THE DEFAULT IS 20 
PERCENT CF THE TOTAL NO OF LINES AND COLUMNS. 

►a******************************************** ******** ************ 


REAL* A F 12 500) .G12500 ),R!2500),ICSETI90) iRNORM.HEANl «MEAN2 
INTEGER**, RUNTA0UO,3 ) , ID 1200) , L DAT A1 ( 2500 ) , LDAT A2( 2500 I 
INI EGER* A LIST1 19) /• *CQR • INPU' , 'BLOC* , ■ F UNC* , * CHAN • , 'SAHP 1 , 

1 'LINE* ,'EXPL" , ' END •/, 

2 LIST2(3)/*RUN • , ' T APfc ' , ' F I LE ' / , 

3 LI ST 3121/' LINE' ,’COLU'/ 

I NT EGER* A ICARO!20),FCH20)» IVEC12) 

INTEGErt*2 CSELISlilCSELUO), ACC I 30) 

LOGICAL* 1 ID AT A 125000) ,FL AG 1 1 Z ) / l Z * . F ALSE. / 

COMMON /CORCOM/ FSTLN ,LL I NE , FSTC L, LCOL »NS I ZEX , NS IZE Y .LACK ,L AGY , 

1 PCTX.PCTY 

FORMAT! IX, 'ERROR IN CTLWRD. ERROR C0DE»',I3,' EXECUTION TERMINAT 

^DRMITIIX, 'ERROR IN CTLPRM. EXECUTION TERMINATED') 

FORMATtlX 'ERROR IN IVAL. EXECUTION TERMINATED') , 

FORMA) (lxI'ERROR IN CHANEL. EXECUTION TERMINATED') 

FORMATtlX, 'MISSING CONTROL CARD OR PARAMETER. EXECUTION TERMINATl 
I'l 


PCTX-20 
P.C T Y « 2 0 
DO 777 1*1,30 
NCC 1 1 l»0 
ICSEL I 1 1 *0 
CONTINUE 

DU 778 1-1,90 
ICSET 1 1 1 —50000.0 

CONTINUE 
LSZ*9 
1ER-0 
INRD* 5 

CALL CTLWKO 1 1 C ARO , 1 C OL ,t 1ST 1 , LSZ , I CODE , I NRD , 1ER I 

GO 1 T0 K 1 9 9 1 101^102?103?IOA.105,106,10T,108»,1 CODE 
FLAG! 1 1 “ • TRUE • 

GO TO 100 

INPUT RUNTABLE 
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CALL 3 CTLPRNUCARO. tCOL, LI ST2. LSI, ICOOE, 10021 

CALL l lVALUCAR0.ICUL,iyEC,LSZ,1003> 

GO TO (501 ,502 , 503 J , (CODE 

run number 

RUNTABll.ll-lVECUl 
FL AGI 2 1 » « TRUE « „„ 

IF t 1C0L. LT.T2I GO TO 101 
GO TO 100 

TAPE NUMBER 

RUNTABt l ,2 I»1VEC( 1 ) 

FLAG! 3I-.TRUE. „ 

IF( IC0L.LI.72) GO TO 101 
GO TO 100 

FILE NUMBER 

RUNTABI l,31*IVEC(i) 

FLAG! A)-. TRUE. „ ... 

I F ( 1C0L.LT.72I GO TO 101 
GO TO 100 


CALL CTLPRMI ICARD, ICOL,LIST3,LSZ,ICODE,IOOZ) 

CALL 2 I VAL( ICAROflCOL, 1VEC ,LSZ,1003) 

GO TO I601,o02), I COO E 

NO OF LINES 

FSTLNMVECU) 

LL I NE * I VEC 12 I 
FL A G ( 5 ) * . TRUE • 

IF I 1C0L.LT.72) GO TO 102 
GO TO 100 

NO OF COLUMNS 

FSTCL* I VEC 111 
LCOL *IVEC(2l 
FLAG(6)».T RUE . 

1F1 ICOL.LT.72) GO TO 102 
GO TO 100 

AUTO OR CROSS 

00 151 I -1.20. , 

FCTII ) * 1 CARD 1 1 I 
CONTINUE . „ 

FLAGI T I «• TRUE. 

GO TO 100 

CHANNELS CARO 

CALL CHANEL l ICARD , ICOLtNCR , ICSEL , I CSET , NCC , 1006 ) 

N01H-NCR 

FLAG! 3)*. TRUE. 

GO TO 100 

CROSS TRACK LAG 

CALL 1 1 VALl ICARO, I COL, I VEC ,LSZ, 1003) 

PCTX-IVECt 1) 

swjMJrtJSVfi G o to id 

GO TO 100 
ALONG TRACK LAG 

^ALL 1 1VALC ICARD, ICOL, IVEC ,LSZ, 1003) 


. . CORELAT FORTRAN B1 

PCTY-IVECI11 
FLAGI 10 1 » . TRUE . 

I F I ICOL.LT.72) GO TO 101 
GO TO 100 

EXPONENTIAL FIT 

FLAGI l 1 ) “.TRUE. 

GO TO 100 

END CARO 
FLAGU2) -.TRUE. 


CHECK IF ALL CONTROL CAROS HAVE BEEN READ 
DO 651 1*1.3 

IFI .NOT. FLAGI II) GO TO 652 
CONTINUE 

I FI .NOT. FLAG 112)) GO TO 652 

GO TO 656 
WRI T E 1 16,10151 
GO TO 999 

CONTINUE 

GO TO 125 

WRI TEI16,101ll 1ER 

GO TO 999 

W-tl TEI 16,1012) 

GO TO 999 
WR! TE 1 16, 1013 ) 

go ra 999 

WRI TEI 16,1016) 

GO TO 999 

NS1ZEX-LC0L-FSTCLH 
NS1 ZEY-LUNE-FSTLNU 


LA3X»PCTX*NSIZEX/100 £3 hh 

LA3Y=PCTY*NSIZEY/100 'N CQ 

CALL CRLT IF,G, R, L DAT *1 , L OAT A2 , RUNT AB, NCC, FCT, FLAG) 


SUBROUTINE CRLT 

PURPOSE PEBF0RM A TH0 DtHENSIL'NAL CORREL A T 1 ON OF AN 1 MAGE IN TWO 
ARBITRARY BANOS ANU DETERMINE A SPATIAL CORRELATION MATRIX. 
DESCRIPTION OF PARAMETERS 


F - ARRAY TO STURF THE INPUT IMAGE (FIRST CHANNEL) 

G - ARRAY TO STORE THE INPUT IMAGE ( SECONO CHANNEL) 

R - RCSULTING SPATIAL CORRELATION MATRIX 

LDA1A1 - STORAGE ARRAY 

L0ATA2 - STUKALE ARRAY 

RUNTABI 1,11 - RUN NUMBER CIF THE UESIRED AREA 

srahiii ■ file ra §p im bribes wii 

per - DE SI^GNAT? NG* AUT 0 OR CROSS FUNCTION 

FLAG - UP T 1 ON ARRAY 

REHARKS 

ARR AY^MUS T "bE^D^ME NS^ONED L I N^IhE^MA?N KSgU’lcg&oKS 

THE T LAG S iN E ESTlMATINi C THE e CORRELATjON FUNCTION IS USER SUPPLIE 

■ SS till IBteWjR&.DE AN 

ESTIMATE WITH SMALL VARIANCE. 



r'ooriorwv onnnnof>onnnonnonr> 
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ARRAYS F.G.AND R ARE EMPTY UPON ENTERING THE SUBROUTINE. 
SUBROUTINE AND FUNCTtON SUBPROGRAMS REQUIRED 
ACF . TAPOP, GETRUN 
METHOD 





SUBROUTINE CRLT IF ,G»R,LDATA1»L0ATA2, RUNT AB.NCC, FCT, FLAG) 
*06606e6«***06**«*» 0 *»**«* *0**06* 6* *0**0* 6**0 666***06*6*0 *****006666*0 
LOCAL VARIABLES DEFINITION 

*500*006** s>**6 0*0 ** 0 <<*1*0 1>* 660**60 60 0*0060 6 60***6066* A **#6***6 66«*#*666 

IMPLICIT INTEGER (A -l) 

RE AL*4 F I NSI ZE X, NS I2EY1.G (NSIZEX. NSIZEY) , R IL AGX. LACY J , 

1. Hfil.RUi ,S(M , MOT 131 ,H« 1 2 1 , WR2 f 2 ) .ERROR (35, 

2 1E«( il.KKSC III ,FN( 31 tFMATI 3 5 .RNORHi ME AN l, PE AN 2, 

^NTEGER'^KUNT ABUuJ 3 )? 1 6^200 > , LOATAl ( NS I ZEX J,L0ATA2IN$f ZEX) 
Integer#* acor/* aut’/.ccor/' cro</ _ .... ..... 

INTEGER*'* FrtTl f 61/ '(///',' 30 X,** 0*,MF5. «. *2 ,IX* M */ 
INTEGER** FHT2C5I/' (',* 0* I FS. », *2,tX* , * >4 •/ 

INTEGER** FCTI2C), Nt /* 0'/ 

I NTEGER62 CSELtet,NCCI30) ... . 

LOGIC AL* 1 10ATAI25G00 1TFLAGI 12 I 

COMMON /CORCOH/ FSTLN,LL1NE »f STCL.LCOLt NS l ZEX, NS IZEY ,LAGX ,L AGY, 

1 ' PCIX.PCTY 

6*#o6ee6e**#66#oo#*6*#6v«*6**oo*66**o*6****##*****6*#**o******66*o6a*# 

COMMON BLOCK VARIABLES DESCRIPTION 

C FSTLN - FIRST LINE 

C FSTCL - FIRST COLUMN 

C LL1NE - LAST LINE 

C LC3L - LAST COLUHN 

C NS l ZEX- NQ Of COLUMNS 

C PCTX -- LAG AS a^pIscentace of nsizex 

C PCTY - LAG AS A PERCENTAGE OF NSIZEY 

£ 0*6*60******** *#» 4 # »»»***# 6 ** 6 *#*#*# 666 * 6 *# 66 6 ** A 0# A 6 6**6 06 *60*66*6**6 

60 FORMAT! I X, *GFTRLN 6RROR**,I2) 

61 F0RMATUX,»10PFS ERRUH«»,lZI 

“ »* > 121 


62 FORMATILR.MOPRV ERROR.', IZJ 

INUNI T«12 
LAG*3 

C* *0*60000* 600 *080** »#0##06 0#6#*#*06***#*0#*#66*0####6***6#***0600*660* 

C* SETUP THE RUNTABLE ANC FIND THE PROPER FLIGHTLINE 

r 

£ 0 ** 60 **** o ##* o #******** eo * 6 e #*# 6 * oo *#» 6 *#** 6 «** o ** o # o » o » 6 * 6 «** ooeooooo 


FSTCL*FSTCL-I 
RUNSEL«RUNTAB( 1,11 


CALL GETRUN IRUNSEL , I NUNI T, I O.ERR.RUNTABr U 
lFfERR.NE.Ol HAITEU6.60) ERR 
lFtERR.NE.OI GO TO B8 t 
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£ NRTRY-2 

£*«*# 0 * e * ****#**# 6 #* ******* 0 *# S 6*6 0 * 666 **##* 06 ** 6 *** 66 06*6 ******* 0 * 6 #*» 

C» POSITION THE TAPE AT THE START OF DESIRED RECORD 

^* 6066 * 6****66 6 * 66*6 66*6 ******6 66 ** 66 * 6066 * *************** ••** 6 ******** 

ICOUNT=F STLN-1 
IF( ICOUNT.EQ.O) GO TO 2U 
CALL TUPFS I INUNI T ■ (COUNT ,ERRI 
lFfERR.N6.0j HR ITE l 1 6 ,61 I ERR 
IF ( ERR.NE.OI GO TO BBl 
211 CONTINUE 

C NSC-IDI6I 

NC -IBIS) 

IFIFCTt3I.EQ.AC0R) DS PL-F ! STCL,.,„, 

IF | FCT 1 3 ! . EQ.CCOR J DSPL«N SC* FSTCL 
DO 11 (.1,8 
11 CSELI I l«0 

C IFIFCTI3). EQ . ACOH ) GO TO 14 

CSELINCC 1111*1 
CSELINCCI2I1-1 
£4 CSELINCC 1 1 1 ) >1 

c* 6 *6 ******** ************* *6*6* *66 6* 6* ******** ***666** ***************** 
C* READ IN THE DESIRED AREA 

£•* 6 ** 6 * * 666 * 6*6 66 6 * * 6 * 6*66 6 *** 66 * * 6 66 ****** 66 * 6 ** 666 * •**♦ 6 ***»*** 6 « 0»6 

c 

DO 30 3=1, NSIZEY 

CALL TOPKV I INDMT ,NSC,ERR, I DATA.NRTRY.NC, CSEL.LNID J 
1FIERR.NE.0I HR ITE 1 16 ,62 ) ERR 
If I ERR.NE.OI GO TO B81 


CALL MOVBYT IIOATA, FSTCL, I .LDATA1.3, A.NS IZEXI 
CALL HOVOYT IIOATA, DSFL.l ,LDATA2 ,3, A,NSI ZEX) 


00. ZO 1*1. NSIZEX 
FI tVJ)*LDATAl I 1 1 
G(I,J)*LDATA2<1> 
CONTINUE 
CONTINUE 


C * ******** S * 6 * * 66 * * 6 * 66 ** 66 ♦* 6 »#*##** 6 * 66 ** 66 *** 6 #** 6 **** 6 #*****#*# 6 #»* 
C* FIND THE MEAN OF THE PICTURE 

£0606***6 *****6* *6 ****** 6*6 »*6« *6* *66 6***#*66# **6**66* 666* 6*66 *6*6*6* 6* 

00 200 J*l, NSIZEY 
DO ZOO [•»! .NSIZCX 
NEANl 'ME ANt *F 1 1 ,41 
ME4NZ »MEANZ*G I i ,J1 
ZOO CONTINUE 

ME AN1 *ME ANl/FLO AT I NS I ZEX* NS I ZEY ) 

ME ANZ .MEANZ/FLOAT INS! ZEX6NS IZEY) 


C * ********* * 6* *6 66 ******* ****** * 6666 ****** * 6 *** 66 **** 6 **** ************* 

£ 6 * 6 * 66 ** ***** *»«**».* 6 * 6 *#***»* 6 * 66 * ***»•*»*« * 6 6 ** 6 * •«•»•»«• 6 •*••****( 

C‘ ’ ' 


DO 255 J* I, NSIZEY 
DO 255 I»I, NSIZEX 
FI j ,J f =F I j ,J 5-MCAN 


FI I ,J )=F I I ,J 5-HCANl 

652 FORMA TUX, I 2 UXU2UX ,2 1 1 PE 11. 4 , IX > I 

255 CONTINUE 
220 CONTINUE 

£ 666 * 6 * 6 6 6 * 6 * 6 * 6 * 66 **** 6666 * 666 * 6 6 666 ****** ******* * 6*6 ******** 66 #*** 0 *« 
I COMPUTE THE AUTO I CROSS) CORRELATION FUNCTION 

C « * 0 * 66 * 6 66 * 60 6 * 0 *» 6 » 6 * 6 *** 6 * 66 006 ** 66 * 66 * 6 ** 60 *** 66 * 6 * 6 #*** 6 **** 6 ** a *> 
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CALL ACF IF,G.R,NSIZEX,NStZEY,LAGX,LAGY) 


£ 

£********♦♦♦*♦♦**♦**♦*♦***************************♦*********************' 

C NORMALIZE THE RESULTING CORRELATION FUNCTION „ . a *’ 

£»********♦**«***♦*******«***************♦************♦*******♦********•; 

RNORM»R( 1,1) 

DO 68 J-l.LAGY 
DO 68 I» 1 ,LAGX 
R(l,J)*Rll,J)/(RNQRM) 

68 CONTINUE 

C IF I .NOT. FLAG (111) GO TO 902 

£********♦****♦*******♦♦♦♦*********************♦*♦**********************, 

C EXPONENTIAL CORRELATION MODEL DEVELOPMENT 

£******* I***************************************************** **********! 

C USE A LINEAR LEAST SQUARE FIT TO THE LOGARITHM OF THE FUNCTION 

IFIRIl,n.LT.R(2.l)l GO TO 901 ; 

IFf Rl 1,1 1 .LT.R! 1,2) > CO TO 901 1 

C 

DO 831 l«l,LAG _ l 

IFIRU.n.LE.OJ GO TO 901, 

FN(II*"AL0GIRI LAG - II-1I.1II 

831 CONTINUE ; 

c CALL LSF IH,WGT,P,S,FN,FHAT,C,IER, ERROR, WRl, WR2 , LAG, RKSE ,RHU> ; 

ERl’HMSEIl) 

RH0X-EXPI-CI2) 1 
C 

00 832 1*1, LAG 


1FIRI 1,1 ) .LE.OI GO TO 901 , 
FNI 1I—ALCGIRI l.LAG-l I— 1 1 1 1 
CONTINUE 


CALL LSF (H, WGT, P, S.FN.FHAT.C.IER, ERROR, WRl, WR2, LAG, RHSE, RMU) 
ER2 *R MSt 1 1 1 
RH0Y*EXP(-C(2) ) 

CD TD 902 
WRITE 116, t0l3> 

H? I J.F 19 l AG 13 I i nncr nccroTDC -rue ATIi no I 


013 FORMATllx, ’EXPONENTIAL FUNCTION DDES 
IF THIS AREA. REQUEST IGNORED’! 


NOT DESCRIBE THE ACF OR CCF 0 


902 CONTINUE 

£*•******♦**♦**♦*♦*♦****♦**♦*♦♦ A************** ******** ****************** 

C GENERATE THE OUTPUT 

£*♦*******«♦♦*,**♦*♦♦*********«•♦************************************♦*** 

1 IF IFC 1 1 3 ) .EQ.ACCR) NCCI 2 ! -NCCI 1 ) 

FMT1I 3 1 * L AGX *N 1 
1FILAGX.GE.10I FHT 1 1 3 !*N1 *9 
FMT2(2I*FHT1I3! 

WRITE 16, IOOO) 

1000 FDRMAT(IhI) 

1001 FORMAT I 40 X?’ TWO DIMENSIONAL SPATIAL CORRELATION ANALYSIS’! 

WRITE 16, 1002! 

1002 NCCU..NCCI 2 ! 

1003 FORMA T ( 50 X,' CHAANELSO ’, 2 ( 12 , 1 X 1 ! 

C (CLASS! I t 2 », 1 * 1 , 2 ! 

1004 F 0 RHATI 40 X, ’CLASS' . 2 * 4 ) 


CORELAT FORTRAN 


WR1 TEI6, 10051 
FORMAT! ////I 
WRl TE (6, 10061 

FORMAT! J0X.’2-0 SPATIAL CORRELATION MATRIX* I 
WRITE!6,1G02) 

00 320 J*1,LAGY 

WR I TE 1 6, FM 1 1 ) IRII.JI ,1*1 ,LAGX) 

CONTINUE 

IF (R( 1,1 ).LT.R(2,1!1 GO TO 671 
IFIR! 1,1 I.LT.R! 1,2)1 GO TO 671 
IF ( .NOT. FLAG 1 1 1 1 1 GO TO 671 

WRITEI6, 10051 
WR 1 TE 1 6, 10 07 1 

FDRMAmoX, ’WEIGHTED LEAST SQUARES FIT INFORMATION’) 

WR [ TE ( 6, 1002 I 
WR I TE I 6, 1006 1 RMU 

F0RMATI30X, ’WEIGHTING MATRIX DIAGONAL BASE*’, FA. 2) 

WRITE !6,‘ 10 02 ) 

WRITE (6, 1009) ER1 

FORMAT! 3UX, 'WEIGHTED LSF ERROR (CROSS TR ACK I * * , E 14. 7 1 

WRITEI6, 10 02 1 

WR I TE ( 6 , 10 10 1 KFOX 

F0RMAU30X,'ADJACENT SAMPLE CORREL'AT ION- • , E14.7 1 
WR I TEib, 1002) 

FORMA itiJXj' WEIGHTED LSF ERROR (ALONG TRACK I *’ , E 14.7 ) 

WR I TE (6, 1002 ) 

WR I TE ( 6 , 10 12 ) RHOY . „ 

F0RHATI30X, 'ADJACENT LINE CORRELATION »*,E14.71 

CONTINUE 
WRITE! 16,1002) 

FORMAT! 1 1 

FORMAT! 15X?°TW0 DIMENSIONAL SPATIAL CORRELATION ANALYSIS’) 
WRITEIL6, 20011 

WR I TE 1 1 6 ,2 003 ) NCC ( 1 ) ,NCC ( 2 ) 

F0RMATI20X ,’CHANNELSO *,2 I 12, IX)) 

WRITE! 16,20011 

F0KMATil}x??2-D SPATIAL CORRELATION MATRIX’) 

WR J TE I 1 6 , 2 001 ) 

DO 330 J * 1 ,L AG Y 

WRl TE ( I6.FMT2 ) I R! I . J ) , I » 1 , L AGX I 

CONTINUE 

WR I TE ( 16 , 1 005 ) 

RETURN 

END 


USAGE 

CALL ACF (F,G,R»NSIZEX,NS1ZEY, L^GX , LAG Y ) 

D p SCR,PT i D ^SY P C a ga? E AlSnG THE AREA TO BE CORRELATE D( 1ST d 
G - ARRAY CONTAINING THE AREA TO BE CORRE LATE D! 2ND Cf 

R - RESULTING SPATIAL CORRELATION MATRIX 

NSIZEX - NUMeER OF COLUMNS IN THE PICTURE 
NS1ZEY - NUMBER OF LINES IN THE PICTURE, 

LAGX - LAG IN PIXELS ALONG THE COLUMNS 

LAGY - LAG IN PIXELS ALONG THE LINES 

REMARKS 


SUBROUTINE AND FUNCTION SUBPROGRAMS REQUIRED 


THE LAGGED PRODUCT SUM METHOD IS CARRIED OUT 


ORIGINAL RAGE 
OF POOR QUALT 
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9 

c ************* 

c 

c SUBROUTINE ACF IF,G,R,NS1ZEX,NSI ZEY, L AGX »L AGY I 
£ LOCAL VARIABLES DEFINITION 

IMPLICIT INTEGER IA-Z) 

1 REM -** ^^^^J^U’ G,flSIIEX * NS '«yl.MLACX.LAGYI f 
ROlH=NSIZExiNSlZEY 
£«»*»* *0 0«"> «*»««»**«*»»« *»«««,« 

c SUN THE LAGGED PRODUCTS 

00 10 £TA*1,LAGY 
DO 10 TAUM,LAGX 
oo 30 j «i , ns tier 
c DO 30 1 oi.NSUEX 

C ,F( * 1 ttTAU-l I .GT.NSUEX1.0R. ( I JaETA-U.GT.NSIZEY )J GOTO Z5 
C . 

25 CONTINUE A ^ ■'RITAU.ETAIaGII +TAU— I » J*ETA-»1 > »F < I , j} 

30 CONTINUE 

io . ?^?y^r , ‘ R,TAu ' ETA,/RD, « 

RETURN 




^ * mw • Q 


SUBROUTINE LSF 


PURPOSE 

TO FIT A LINEAR FUNCTION THROUGH A DATA SET.USING LEAST SQUARES 
DESCRIPTION OF PARAMETERS 

m^?,s f g mn 0,5 UNEAft FuNcno « s 

P -HTRAN$POSE*HGT 

R -P*H 
S -RINVERSE»P 
C -LEAST SQUARES COEFF. 

U^KHriW* VAt - U£S 

ERR-HE IGH TED F-FHAT 
LAU-Nl) OF DATA VALUES 
WR1-WURK VECTOR 
MR2-MQrtK VECTOR 


REMARKS 


THE EQUATION CF THE LINEAR FIT IS OF THE FOKN.fl.rj** 

the following HATRich are oimInsioneo ZalAG* 1 2 x 

H.P.S 

THE FOLLOWING HATfilCE S ARE OIHENS IONEO LAf IN THE MAIN PROS 
UGT.HRl »WR2 .ERROR, ERR, FN, FHAT 
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THE FOLLOWING MATRICES HAVE CONSTANT DIMENSIONS 
R!A> ,RKSEI IlfCIZI , 

SUBROUTINE AND FUNCTION SUBPROGRAMS REQUIRED 
TPRO.GHPRD.HINV 


METHOD 


?^I 5 ^l 23 V^fffI^SSIM eCHUIOUe > S I^EMCNTFO.FORMULA USE 
F ,°5 E^K2F R .S£I AILS SEE 1 DICRETE PARAMETER ESTIMATION *,6Y ( 

J MtfiOtL « { 973 ( 


SUBROUTINE LSF IH, HOT ,P,S ,F» FHAT, C, ERR, ERROR ,WRI ,WRZ , LAG. RHSE ,RHU) 

*0 «* £64*6* *Crt D*«i5.**UiiJ,<» <«*♦♦»»» 

C LOCAL VARIABLES DEFINITION 

C* ****«« ********** «« «»«*««* »**o ♦„**.; »**.<,»» 

*’ IMPLICIT INTEGER IA-ZI 

c t R6AL * 4 ,u4ni??smi?uA5iiH:^uur45:^i* Cfi ** eRRii, - eftRoaci > 

80 FORMAT!//! 

81 FORMAT!///) 

NPNTS-LAG 

DO 30 1 - 1 iNPNTS 
H!I (»1 

30 CON T I HUE 

K1 « NPNTS ♦ I 
X2»Z*NPNTS 
DO AO I * K1 ,KZ 
H! I)»KZ-I 
AO CONTINUE 

C*p«aa»*ev««**«*e»*a*»**<t»« .(it.ti,.,,,,,,,,.,,, » ***°«e a *o* *»*.***» **•«««• 

C» INITIALIZE THE WEIGHT ING MATRIX 

C*o** »* P **<.«*» o*se e» •« *» ««*** «»« »»*** .... ******** ** ***** ******* 


, RHU-O.A 

C 

00 28 I >1, NPNTS 
K*NP.NTS-1 
WGTIt )=R«U**X 
28 CONTINUE 

C 

£«***«*******••******•«*•****** 

C DETERMINE THE LEAST SQUARES FIT COEFFICIENTS 

CPP****** ***** ************* **** ****»«*•**#***« ******** *********** ** ***** 

c 

N-NPNTS 

M«Z 

MSA-0 

MSB-2 

L-NPNTS 

CALL TPRD IH.WGT,P,N,M,MSA,MS6,U 

V 

H-2 

H-NPNYS 

*-■2 


•w. 1 



FILE. 


. . CORELAT FORTRAN B1 
CALL GHPRD ( P.H. R. N.H tL) 

CALL MtNV (R,N.D£r,HRl,HR2> 


N*2 

H »2 

L»NPNTS 

CALL GHPRO IR,P.S,N,H,U 
N *2 

H»NPNTS 

CALL GHPRD t S.F »C»N» H ,L> 


CD) 

COl 

CO< 

CUI 

COl 

COl 

cm 

cui 

cm 

cm 

coi 

COl 

cm 
co i 
coi 
cm 


r******************** ****** **** ********** ********************* **********COI 

£I.***«J!!!*«™!;A*2******J*2***U**I*******************«*****************COf 


N'NPNTS 

H-2 

CALL GHPRD CH.C.FHAT.N.K.L) 


cm 
coi 
COI 
COI 
cm 

------- cui 

y coi 

c, ,*******♦»***♦♦**♦♦**♦ a*****************#*****************************^! 

^:,*.*5i?5.WI.S!iiSS^S.k 6 .^ s U52ili55l*^525 ********.************ftj 

' DO 100 1 ■ 1 .NPNTS SR! 

ERRORUI»FUI-FHATin 5UI 

100 CONTINUE 


N-NPNTS 
H« 1 
H5A»0 
H$B»2 

CALL N TPR0 (ERROR. HOT i ERR.N.M.HSA.MSB.D 
N«1 

M«NPNTS 

CALL GHPRD I ERR .ERROR .RHSE.N.M ,L 1 


RETURN 

ENO 


CUI 

cm 

cm 

cm 

CUI 

cm 

cm 

COI 

COI 

COI 

COI 

cm 

COl 

CDI 

COI 

COI 

COI 

COl 


ORIGINAL _L 
of poor quality 
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C SUPERVISOR FOR THE COMPUTATION OF THE SCANNER OUTPUT STATISTICS 

C WRITTEN 09/25/78 BIJAN G. KOBASSER! 

c 

<5*o»ooo4*c4o4o<i»*404*«44o *44»ti**44*4sea*444*44**»#«i444*»a**4*a**44t!i»ij* 

OESCRIPT ION ANO PURPOSE 

THE SCANNER OUTPUT STATISTICS PROGRAM IMPLEMENTS A LINEAR 
TRANSFORMATION TO OBTAIN A SET OF STATISTICS AT THE OUIPUT 
OF A MULflOANO SCANNING RADIOMETER IN TERMS OF THE CORRESPONOIN 
INPUT QUANT HIES. IHIS TRANSFORMATION IS ACCOMPLISHED BY 
THE COMPU f A I ION OF THE * SCANNER CHARACTERISTIC FUNCTION AN 
ANALYTICAL EXPRESSION. THE MAIN PARAMETERS OF WHICH ARE THE 
IFOV SUE AND INFORMATION ON UATA SPATIAL CORRELATION, 

DESCRIPTION OF CONTROL CARDS 

RSC6NSTAT 

THIS CARD SPECIFIES THF PARTICULAR PROCESSOR REQUESTED 
CHANNELS 

fHE DESIRED SUBSET OF THE AVILA81E CHANNELS IS GIVEN HERE. 

IT IS IMPORTANT TO REMEMBER THAT THE NUMBERS APPEARING 
ON THIS CARO ARE THE OROER OF THE SELECTEO CHANNELS NCT 
TMthR ACTUAL NUMBER. FOR EXAMPLE IF THE AVAILABLE CHANNELS ‘ 

ARE 8,9,12*14 ANO CHANNELS 8,9 ANO 14 ARE SELECTED, CHANNELS 
CARD SHOULD READ t ,2, A 

CLASSES 

THIS CARO SPECIFIES THE NAME OF EACH CLASS. EACH NAME MUST -BE 
PLACED IN A FIELD 7 CHARACTER LONG FOLLOWED BY A BLANK. THE 
CONTINUATION CARD. IF REQUIRED, MUST HAVE THE WORD 'CLASSES* 

IN THE BEGINNING FOLLOWED BY THE' REST OF THE NAMES. 


THIS CARD SPECIFIES THE SPATIAL RESOLUTION OF THE OUTPUT DATA 
IN TERMS OF THE INPUT. BASICALLY IT IS THE NUMBER OF HIGH 
RESOLUTION PIXELS HITHIN ONE IFOV OF THE SCANNER. E. G. FOR 
MSS OPERATING ON A 6 MFTER DATA AND IF0V*5 ACCORDING TO THE 
ABOVE CONVENTION THE ACTUAL SPATIAL RESOLUTION IS 30 METERS 

APERTURE 

THC CHOICES HERE ARE 'GAUSSIAN* OR 'RECTANGULAR* 

SNR (OPTIONAL I 

THIS- CARD SIMULATES THE EFFECT OF RANDOM ADDICTIVE NOISE ON THE 
POPULATIONS SI ATI S TICS AT THE SCANNER OUTPUT.THE NOISE COVARIAN 
MATRIX is DIAGONAL WITH OFF DIAGONAL ELEMENTS EQUAL TO ZERO. 

THE SNR IS DEFINED AS THE RATIO OF SIGNAL ENERGY (DIAGONAL ELEH 
ENTS) TO NOISE ENERGY. SNR MUST BE GIVEN IN DECIBELS DEFINED 
SNR » IQ.*ALCOtO( SIGNAL ENERGY/NOISE ENERGY). THIS SNR WILL 
BE THE SAME IN ALL CHANNELS. 

PUNCH IDPTIONALI 

IF PRESENT THE OUTPUT STATISTCS IS PUNCHED OUT. NOTE THAT SOME 
BLANK CAROS ARE INCLUDED tN THE DECK FOR COMPATIBILITY REASONS 


THIS ARD. SIGNALS THE END OF CONTROL CARDS. DATA FOLLOWS IMMED1A 
TELY 

INPUT DATA STRUCTURE 

INPUT TO *SU$ CONSISTS OF 3 SEPERATE DECKSO , , „ „ 

i- LARSYS STATISTICS DECK WITHNCJ CHANGES IT SHOULD 
HOWEVER Bf. In 44 CHARACTER FORMAT ** 
l- SPATIAL CORRELATION PARAMETERS ARE ENTERED VIA> A NDIH X' NDIM 


FILE. . . SCANSTAT FORTRAN 81 

C MATRIX. THE 11,41 ElEHENT OF IT IS THE PIXEL-TC PtXbl 

\ ^E RR E{iSU D ?R^K^^L L Ir{or IH,S DECK “^SPONDS »« 

C 3- SAME AS 2 EXCEPT FOR ALONG TRACK DIRECT ION. 

C REMARKS 

9 IHA s » p 5S52^,. ls „5 u S5 E>i ' ri - y capable of processing up to zo classe 

i g?,$, sto sm r i8 »*sm r «8 M!i&,??ii s )a? s 

l i!8r'85IS5%‘S,55?abI",< u 51!I» ! !5 , If,)'5I-SS c EoESSSK tD 

C HOW TO RUN THE PROGRAM 

C THE SOURCE ANO TEXT F JLES ARE LOCATED ON THE DHDSK ANO OHSYS 

C DISKS RESPECTIVELY. AUTHCRIZED ID'S ARE AUTOMATICALLY LINKED 

C TO BOTH DISKS AT LOGIN TIME. OTHERWISE ANYONE CAN ACCESS THE 

C DISKS THRU IKE COMMANDS 'GETOlSK DHSYS' FOR TEXT AND 

C 'GtTDISK DHDSK' F CR THE SOURCE, 

C AFTER THE PRUPER LINKS ARE ESTABLISHED, TYPE IN »SOS. 

C EXAMPLE OF THE CONTROL CARO SET UP 

C 4SCANSTAT 

C CHANNELS l, 2, A 

C CLASSES BARE SOI CORN SOYBEAN WHEAT 

C IFOV 2 

C APERTURE GAUSSIAN 

C PUNCH 

C , END 

, 1ST EGER* A LIST ( 8| /•« SCAN 'CHAN',' CL AS' I FOV, • A PERN * SNR *, *PUNCN 
1 'END */, I VEC (l).iCAKD120) . NACRI 7201 

INTEGERS BLANK/* V.F STCRD/ 'LARS’/ 

INTEGER* A SD I M , S1GMAX , S GMAY , HE ADI 20 I , APERT t 20) 

INI EGER*A FHT 1(5)/' IT*,',’ ON NFS. N« 2, IX N ’ )l« /, 

1 FHT2 I6J/N75X', NN' ON ' IF8. N » 2, IX ' I I */, 

2 FHT3I6)/*It •,*;* ONMFT.N'ljIxNN ) $ / 

3 FMTA (<.)/* (ACX*. 


FHTAI6)/'(ACXNNN* 
NL/' 0'/ 


1NTEGER*2 SPREAD! 30) , ICSE 1 110 1 . NCI30 ) .NCC I 30 1 
REAL*A CUV IN 1 7 2C ! • COV CUT I 720 ) . SUBCOV I 720 ! , SLPXI ? 20) , SLPY1 720) , 

1 SUBSLX(720jjSUBSLYt7HC),PXj720(iPYl720)|sPXI7 20) , 

2 SPYJ720) ,Cl l720),C2f 720) ,C 720) ,RHOj 7201 , 


£ arTl f £\j I *u i 1 i c\j t t * £ \ f 

3 MU(lbO) i SUUMUI I&OI 

RF AL*4 ICSET 1 90) 

LOGICAL*! FLA0I8S/B*. FALSE./ 


COMMON /SOSCOH/ NC LS.SDIM, NDIM, NS IZE, (SIZE, NTS, NSA,S IGMAX ,S IGHAY , 
1 NT H» N S Nrt 


ICIO FOR HATH x, 'ERROR VN C TLWRD, EXECUTION TERMINATED.*) 

1011 FORMATITX, 'ERROR IN IVAL. EXECUTION T FRM I NAT ED. • I 

1012 FORMAT! IX, 'ERROR IN CHANbl. EXECUTION TERMINATED.') 

1013 FORMATITX, 'MISSING CONTROL WORD. EXECUTION TERMINATED.') 

1015 FORMAT !:IH l ) 

1 0 16 FDRH4r:i7T 

1.017 FORMAT 1 35X, ' S CANNER OUTPUT STATISTIC S») 
1G1’8 F0RMATI20X, 'SCANNER OUIPUT SPECTRAL STATISTICS* 1 

1019 FORMAT i 3SX , * ARER.IuRf 0 3 AA ) 

1020 FORMA T { 2 OX APERTURE 0 * ,3 AA ) 

1021 FORMATIJSX.'IFOV SIZEC',12,* HIGH RESOLUTION PIXELS*) 

1022 FOR MA T ( 2 OX,, ' IFOV SIZEO',12,' HIGH RESOLUTION PIXELS') 

1023 FORMAT! A9X, 'CLASS N 2 AA I 
1Q2A FOR MATI25X, 'CLASS N2AA 

1025 FOIMAmOX,' INPUT COVARIANCE HATR IXNA5X ,* OUTPUT COVARIANCE MATRIX 

1026 |FORMATUX, 'INPUT COVARIANCE MATRIX* , 10X, ' OUTPUT COVARIANCE MATRIX* 
£02 7 FORMAT!///) 

£ DECODE CONTROL CARDS 
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»•,)»»»•«, a,*,,,,,,*,* ««*•♦**♦♦* •••••®*° # **®®****®®**"®**®******** # ®“* 


SUQRQUTI NE SCANER 
PURPOSE 

TO COMPUTE THE SCANNER CHARACTERISTICS FUNCTION AMO GENERATE 
THE TRANSFORMED STATISTICS 

DESCRIPTION OP PARAMETERS 


APERT - ARRAY CONTAINING THEIFOV, SHAPE 

Cl - ONE D I PENS ION AC SCANNER CHARACTERISTIC FUNCTION 

C2 - ONE DIMENSIONAL SCANNER CHARACTERISTIC FUNCTION 

C * - Cl*CZ 

COVIN - INPUT COVARIANCE MATRICES 

inr : flag U 1n C conircl c cardI ,C£S 

NAOR AOOKEsf^OF^HpHATRIX ELEMENTS TO BE DELETED 

NC - THE CHANNEL SET THAT IS NOT REQUESTED 

NCC 7 COMPLEMENT OF NC _ , 

RHCi - SPECTRAL CORRELATION HATRJK 

P* - AUX ARRAY 

SLPX - ^PAHAL A CORRELATION PARAMETERS, CROSS TRACK 

SLPV - SPATIAL CORRELATION PARAMETERS, ALONG TRACK 

SUBSLK - SUflSET OF SLPX 

luucuv - sueslr of the input spectral cov matrices 

SUBHU - SUBSET OF THE INPUT MEAN VECTORS 


REMARKS 


SUBROUTINE AND FUNCTION SUBPROGRAMS REQUIRED 


THE REQUIRED SUBSET OF THE INPUT SPECTRAL AND> SPATIAL’ *CORRELA 

ha! Rices'll com? lteo , . t,hc s canner character ist i c. u nlt.ion is 


Kio T MIB T i!F TA lia 54 i!E§?? fflt P 55 Hhl&%Sf^c| 4 i| 3 SjlKTi» 

gUSI'SSSftMRS THe'eQUALI I Y^CF^INPUT ^'ANO*GUPUT^ ME All ^VECTORS 


SUBROUTINE SCANER 


ICO VINi COVCUT, SUBCOViSLPK, SLPY ,SUBSLK, SUBSLY • PX, 
PY ,C l IcZ ,CinJ t NCC’, RMOfNADR, FL AGtHUjSUBMUg APERT I 


LOCAL VARIABLES DEFINITION 


REAL-H COV ININSA) tMUI'NSMIl if OVOUTINTSI sSUOCOV {NTS J , SUBHU INTO JpSLPKI 
I NSAI.SLPYINSA) >SUB 5LX INTSI » SUB SLY I NTS J (Pit {NTS IVPVINTSl » 

Z ClINnijCZiNTSi.cTNTShRHOlNTi} 

INTEGER? A £811 foil! <£»•! r-! D| J *!|*J * J *|| •{ 


INTEGER*? FMT2 1 A 1/ ' 1 2X, • » - nil - 

WxlMUt SUNA« t §l£A^SOII>^SfisA>.^ftTI»» 

ami nKSH. * 

LOGICAL*! FLAGS O'J 
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o-x 
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snNiiNoi 


ji+sgH}Adis*(x mssns 
, U*SQWjxd1S»{xixiS8nS 
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30N11N03 
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£••#•********* **************************** 4 ******•**•**••« «•***»•****«* 

C 

DO AS 6 JCLS-I.NCLS 

c hds-ijcls-*ii*isi2£ 

00 A55 K01N“I » NOI H 

NAU«K 0 lH*(KDIM*n /2 

VARNSE*< I 1C)**I-SNR/ 10 . } ) *COVOUT 1 HOS+NAD) 
COVUUTtMOS*NAD>«COVQUT1f'OS*NAD)*VARNSE 

*55 COST I NOE 
A56 CONTINUE 

IF I .NOT • FLAG (7 1 I GO TO 5A A 

c 

£*«*♦**♦* ***** ******************** ************ ******** ft**************** 

| PUNCH OUT THE SCANNER OUTPUT STATISTICS 

£********•**#******♦«********** ************************ ft*************** 

c WRI TE(7>508t 

HA I TE ( 7, SCSI 
WRITE ( 7, 523) 

S2J FORMAT IIX.* *» 

WR l TE (7, 510) NCL5,NFL0,NDIH 


NUN 2 » SO I HA l 
03 53 3' 1*1 t'NUNZ 
WRITE (T»523l ’ 
CONTINUE 

FHm2)*H3»N0fH 

FHT2i2l»N2 

KRtTEt7,FMTl) SUBHU 
HR I TE 4 7i FHT2 ) CjCVOUI 

CONTINUE 

RETURN 

END 


SUOROUUNE ACRES 
PURPOSE 

FIND THE DESIRED SUBSET OF A COVARIANCE MATRIX 

DESCRIPTION OF PARAMETERS 

NADR - ADORE SS ARRAY OF THE ENTRIES TO BE DELETED 
NC - ARRAY OF THE CHANNELS TO BE DELETED 
SOIR - TOTAL HO OF BANOS AVAILABLE 
HO !H - DESIRED SUBSET OF SPIN 


REMARKS 

NONE 

SUBROUTINE AND FUNCTION SUBPROGRAMS REQUIRED 


METHOD 


USING A ONE DIMENSIONAL STORAGE MODE FOR EACH COVARIANCE 
MATRIX IN AN UPPER TRIANGULAR FORM, THE PARAMETRIC ADDRESS! 
FOR EACH ENTRY IS DERIVED. USING THE NC ARRAY, THE LOCATIOA 
UF ALL THE ENTRIES THAT LIE IN THE UNWANTED LINES AND CQLUf 
ts COMPUTED AND STORED IN NADR. 


nnno 
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SCANSTAT FORTRAN Bl 


C 

C 

c 

c 

c 

c 

c 


SUBROUTINE ADRE5 INAOft,NC,SDIM,NDIM) 
INTEGER#* NADR H } ,PD I M,$0 IN 

INTEGERS NCUI 

PD I M» SOIM-NDIH 
lF-4PDIM.tE.0i GO TO 100 
DD 20 J* I • POI H 
Nl*NCl J) 

00 10 1*1, Nl 

NADR I KT* 4 NC t J )* I NC4 J l-I 3/2 * ♦ I 


CONTINUE 


CONTINUE 


DO *5 J-l.PDIH 

NFlN»SOIM~NC< J)H 
DO 35 I » ItNF IN 
K*K* 1 

NADR IK)" 1 1 MC( J )*! “i )* (NC IJI + 1-2 )/2 )#NC( J I 


CONTINUE 

CONTINUE 


RETURN 

END 


fO 



5PRDCT FORTRAN PI 


SPRDCT 

PURPOSE 

EXOSYS DATA IN PUNCHEO FORMAT IS READ AN3 STORED ON TAPE 

REVISED 

3 JULY, 1978 


COMMON tOUOOJ 

INTEGER* 4 IN(5>iOAY<3lfTIM6l3) ,NOGRPS .GROUPS ( 5 ) 
INTEGER*^ I SAM. SINT, D ATE ( 151, INFO! 10,7) 


nicnoiM ijnn.su, i , u h i r i m.inruuut . 

REAL*4 WUC0£F<2,5) , DA TA 1 2 500 1 , X < 100) 

E3UI VALENCE! DATE II i, IO( If ),UNF0|1, 11, ID! 30) ),UNF0! 1,2) , 10(401), 
*HNFO( 1.3), 10! 50ii, (INFO! 1, 4 ), 10160) M l NFO! 1, 5 ) , ID! 70) > , UNFOl 1 , 
• I. 10180)1, HNFOU, 7), 10(901) 

REWIND 11 
NT * 100 

ID INFORMATION 

MR ITE (16, 10) 

FORMA T(5X, ‘TYPE IN DATE**/ IX ,15 (*/•)) 

READ! 15, 15 IDA TE 
FORMAT I i 5A 1 ) 

EXP. NO., NUMBER DFCL ASSES, AND NUMBER OF DIMENSIONS 
WRITE! 16,20) 

FORMA T!5X,' TYPE EXP. NO. .CLASSES , AND DIMENSIONS',/' /// U III' 
READ! 15,25110! 16), ID! 1 7), ID! 18) 

FORMAT! I 3»2X» 12 ,3X, I 3 ) 

EXPERIMENT INFORMATION 

NCLS * 10(17! 

DO 35 1*1, NCLS 
HRITE(16, 10) 1 

FDRMAT(5X,‘TYPE CLASS INFO AND NO . SAMPLES FOR CLASS *, It, /IX, 10 
AX, •///•) 

READ! 15, AC) I JHFOIl , I),L *1,10), 10(20*1) 

FDRMAT(I0Al,4X,13) 

MR I IE ( 1 1 i 10 
CALL SPL8L 

READ SAMPLE FUNCTIONS FROM EACH CLASS 

DO 500 K»1 ,NCLS 
WRlTLl6,aoi 

F0RHAT15!/ I, IOX, ‘SAMPLE FUNCTIONS*) 

NF » 10 1 20 'K ) 

00 ZOO JJ=1.NF 

0 REA0C5,1OOQiOAY.TIME, JNjNQGRPS jNOSAMS 

00 FORMAT? JA4,2X, 3 (12,1X1, /, ASX. 5AA,/,20X,| l, 8X, 13, /> 

RE AO (5, 1 1001 ( DATA! t),I»l, NOSAMS) 

DO FORMAT i 20A4) 

WRITEI6,150)IN 
0 FORMAT! I0X.20AA) 

03 160 1*1 «NT 
3 Xi IT « DATi(I*ll 
MR ITE! 11) X 
0 CONTINUE 
0 CONTINUE 

stdp f1l£ “ 

END 


1 KID! 30) )» ! I NFO! 1,2 > , ID! AO) ) , 
) i , ( 1 NFO! l,5),ID!70)>,UNF0tl, 


non co noovft onoononnnnnnnnnnnoonnnononnnnonn 
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SPQPTM FORTRAN PI 


SPOPTM 

PURPOSE 

TO DESIGN THE OPTIMUM SENSOR FOR A GIVEN DATA SET, 

USAGE 

CALLED FROM EXEC ROUTINE 

DESCRIPTION OF PARAMETERS 

AM - MEAN VECTOR OF DATA 

COV - COVARIANCE MATRIX OF DATA 
PHI - MATRIX OF EIGENVECTORS. 

GAM * EIGENVALUES 
N - DIMENSIONALITY OF DATA SET 
NCLS - NUMBER OF CLASSES 

SUBROUTINE AND FUNCTION SUBPROGRAMS CALLED 
EIGENP.EISORT , 5PWGT1 

METHOD 

THE KARHUNEN-LOEVE EXPANSION WITH THE MAXIMUM LIKELIHOOD ESTIMATE 
OF THE COVARIANCE MATRIX AS THE KERNEL IS USED TO REPRESENT THE 
RANDOM PROCESS. 

REVISED 
14 AUG, 1976 


common itmooi 

REAL* A AH! 100 ), Y(IOO) , COV < 5050 ), PHI P ! 100, 100 ) 

REAL*8 VEC 1 ! 100, 100) • tVI! 100 ), INDICT 100) , ACOVC 100,100) * GAM! 1QQ) 
RE A L* 8 PHI (100,100), SUM 
REAL* A XllOOl , H( 100) 

REWIND 2 
MR i (E( 16,51 

FDRHAT(5X, 'OPTIMUM SENSOR DESIGN') 

READ ID INFORMATION 

READ(2) 10 
WR ITE(6, 8) 

„ F0RMAT(Ihi,5(/)I 

call spldl 

N * IDIIB) 

NCT • N* 1 M* l ) / 2 

COMPUTE COVARIANCE 


WRITE (16, 10) 
FORMA T t SX , 'CO 
NCLS « loll?) 


COVARIANCE BEING ESTIMATED (SPOPYM)') 


DO 20 I* I, NCLS 
NFT * IF T * 10120*1) 

CON * OFLOATtNm/OFLOATlNFT-t) 

DO 30 I* L ,N 
AMID * 0,0 
00 35 1*1, NCT 
COV! I i * 0.0 
00 65 IJ-UNFT 
RE AD t 2 1 X 
IN » 0 
DO 50 1*1, N 

AMU! * AHtll ♦ X! I )/OFLQAT (NFT) 

DO 50 J*l,l 
IN * IN ♦ i 

COV (IN) • COV(IN) * X I ( )P X 14 l/OFLOAT (NFT-1 ) 

CONTINUE 

CONTINUE 

IN • 0 

DO to 1*1 , H 

00 60 J«l, I 

IN " IN ♦ 1 

COV I INI • COVIINI - CON* AM 1 1 )*AHt J) 


CON* AM 1 1 )*AHt J) 


i - 1 

Ni 
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>0 CONTINUE 
I WEIGHTING FUNCTION 
IN - 0 

DO 210 I«l,N 

DO 210 J-l.I 

IN - IN ♦ I 

ACOVt ! - CDV I IN > 

ACOVt J, I i ■ COVUNI 
!10 CONTINUE 

[ CALL SPWGT3IW) 

00 250 1*1, N 

DO 250 J«1,N . , , 

250 ACOVt 1 1 J I « ACOVII ,J)*Wt J) 

■ COHPUTE TRACE OF COVARIANCE 

SUM *0.0 
DO 80 I * 1 1 N 

90 SUM » SUM ♦ ACOVU.I ) 

" COMPUTE EIGENVALUES AND EIGENVECTORS 

75 FORMA T ( 5X i 'E I GENVALUE S AND EIGENVECTORS t E IGENP ) * > 

NM - N 

CALL 5 EiGENPlN,NH,ACOV,TrGAM,EVIiPHI»VECI.INDIC,HI 
CALL E i SORT t N, GAM t PHI ) 

C 

C PRINT EIGENVALUES AND MEAN-SQUARE ERROR 

CO » FLOAT INFTI/(FL0ATtNFT 7 llJFL0ATINFT-l>| 

Cl » FLOA T t A*NF I— 1 1 / ( FLOAT (NFT-1 1 *FLOAT( NFT-1M 

110 FORKAnU/llsX^N* . 5 X , * E I GENVALUE' * 5X i 1 V AR I G AN) • *5X, 'VARIPHI)* »5X 
••MEAN-SQUARE ERROR 1 ) 

OD 150 1-1,30 
VARP « 0.0. 

DO 120 J-1,100. 

VARP - E VARp' t C C0*GA?1l l)*GANt J)/(GAM,tl) - GAM(J),1**2 
115 CONTINUE' 

120 CONTINUE 

VARG « C-l*,GA,Mt I l*GAH ( I ) 

SUM - SUM - GAMUT 

WRITE 16, 1A5)I »GAM( 1 ) • VARG , VARP, SUM 4 ,, 

IAS FORMA Tl AX, I2iAX,F10. A. AX, FLO .A , 2X, F 1 0 . A , 2X , F 1A . 6 ) 

150 CONTINUE 

00* 155 J«t,N 

00 155 l-l.N, . 

1,55 P.HIPt l,J| » PHM*»JI 
DO 180 J»1 ,20- 

WRITE 17,1601 (PH IP 11,0), I -1,N) 

160, FORMA T 1 20AA) 

180 CONTINUE ' 

STOP. 

END' 
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FORTRAN pi 
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COMMON 10(1001 

SELECT NUMBER OF TERMS 
HR I TE ( 16 , 101 

£issfi5;n;«vgj}^ 0F terms •» 

FORMA T ( I 2 ) 

REMIND 2 
RE AO! 21 |0 
WCLS - 10(17) 

N - IDT 18 ! 

NOT * N TERN* IN TERM * l>/2 

NF T * 0 

00 20 I* l, NCLS 

86 T » 

c^U,^- /aaA ' ,mCLS » 

FORMAT(!oFllA) , ^ tl * I *^ ,Ct ' S ^ 
COMPUTE MEAN FUNCTION 
„ 00 300 I*l,N 

0 AMI 1 1 » o.o 

00 320 K-l.NFT 
RE AO 1 2 ) K 
00 320 1*1, N 

0 CONTINUE 4,1 *’ * X «‘^' : ‘-«T(NFT) 

REMIND 2 
RfAO(Z) ID 

READ EIGENVECTORS 
00 40 J-l.NTEBK, 


lkj *g jii.niCKn 

READ! S.Jsl (PHI 1 I, J),| *l,N) 

F0RM4 T | 20A4) 

CONTINUE 

LOOP ON THE SAMPLE FUNCTIONS IN THE DATA SET 
DO E |oo"lCLS=l,NCLS 

M.ici§n R 2.o 
UyVi.l'ztiTl o.o 

NF - I0I2C4ICLS) 

CON ; FLOAT (NFI/FLOAT INF- l ) 

03 ISO I SAM-1, NF 

RfcAD SAMPLE POINTS FROM FUNCTION 
READI2) X 

TRANSFORM DATA USING BASIS FUNCTIONS 
DO 70 J»l,NTERM 

vi j i • o.o 

PQ 70 I- I ,N 

TIJI - YlJI ♦ PHI (1, 3 IP 1X11) - AMI It) 
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